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DIGITAL PROCESSING OF RADIOGRAPHIC IMAGES . 


SUMMARY 


In the nondestructive testing of metallic structures using radio- 
graphs it is sometimes necessary to enhance the images to aid visual 
interpretation. The digital computer is a very useful tool for this 
purpose, providing a means to the accurate and flexible implementation 
of well known filtering techniques and to the development of new tech- 
niques . 


This report provides a description of some of the techniques for 
image handling and image processing on the computer. Both the frequency 
and spatial domain designs and implementations of digital filters are 
covered. Examples of digitally processed synthetic test patterns as well 
as radiographs are shown to illustrate the techniques. The software 
developed implementing the processing algorithms is documented. 

It is seen that digital processing enhances the radiographic images 
quite noticably, but noise in the radiographs is a limiting factor. 

Further work needs to be done in the application of feature-oriented image 
enhancement and pattern recognition methods to radiographs. 
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SECTION I 


INTRODUCTION 


A. Radiography and the Digital Computer 

The basic problem considered in this report is that of facilitating 
the identification of internal flaws in metals by using digital computer 
techniques of image enhancement on the radiographs of these metals. 

Conventionally, the hidden flaws in metals such as blow holes in 
a welded joint and hairline cracks are detected visually by looking at 
dense radiographs of the suspected areas under a bright light with a 
magnifying glass. While many of the flaws can be identified in this 
manner, some cases arise in which the existence of flaws determined 
visually becomes questionable. In such cases it is necessary to enhance 
the radiographs by some method to confirm the existence of flaws. A 
typical example of such a problem is one in which the beginning of a 
crack is visible in the radiograph, but the end of the crack is very thin 
and merges into the background of the radiograph so that it is difficult 
to find how far the crack really extends. 

A large variety of techniques exists for image enhancement. These 
techniques can be grouped broadly into two classes, viz., optical and 
digital. Optical methods consist of preparing transpariences which 
represent desired filtering operations and producing images by passing 
coherent light (or non-coherent light in some techniques) through the 
given picture and the filter transparancies . These methods handle all 
points in a given picture at one time and hence are virtually instanta- 
neous. They are, however, restricted mainly to linear operations and 
their accuracy is limited by the accuracy with which physical filters 
(transparencies or apertures) can be made. On the other hand, digital 
techniques are very flexible and accurate. They can be used for nonlinear 
operations also in a controlled manner. In fact some operations such as 
adaptive filtering where the filter characteristics change according to 
the data content of the part of the image being processed are impossible 
using optical methods exclusively, are extremely cumbersome using electro- 
optic methods, but can be performed digitally with comparative ease. 

The main disadvantage of digital methods is that they handle the picture 
data sequentially and therefore they are slow. But, their accuracy and 
flexibility render the digital computer a very useful tool at least 
during the technique development stage when one wants to explore various 
processes of image enhancement. Many of the techniques so developed can 
in fact be implemented optically ; for fast processing on a production basis. 

We shall concentrate here on the digital techniques since at the 
present stage of learning and developing enhancement methods for radio- 
graphic images it is necessary to use digital methods. 
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B. Definition of Digital Processing of Pictures 


1. Picture Function . A "picture" can be mathematically describ- 
ed by a real function S (x, y) of two real variables x and y which gives 
the "gray level" of the picture at the point (x, y) « Figure 1„1 shows a 
picture whose picture function is given by 

s (x, y) = 1 for (0<x<_.2, 0<y^..5) 

(x= .5, 0<y<„5) 

(•8<x^l, O^.y^.5) 

(0!x<l, .5<y<l) 

and 

S (x, y) = 0 elsewhere 

2. Noise and Filtering . In the process of producing a picture 
there are several factors besides the true image itself that contribute 
to the final picture function. Some of the factors are [l] 

(a) Random noise in the film 

(b) Background of the object surrounding the desired 
part of the object 

(c) Fluctuations in the light output of the image 
scanning system 

(d) Electrical noise in the output of the light 
sensing device in the scanner 

(e) Finite aperture size of the scanner 

All the contributions to the picture function other than due to the true 
image (or the desirable part of the true image) will be termed noise. 

If noise is additive and the "picture function" of noise is v (x, y) , then 
the actual picture function obtained will be 

s ! ( x > y) = S (x, y) + v(x, y) (1.1) 

When the desired picture function is S (x, y) , in general the actual 
picture function obtained can be written as 

s ! (x, y) = f(S (x, y), v (x; y)) (1.2) 

The process of improving the quality of the image consists of filtering 
out the noise v(x, y) so that S (x, y) can be recovered from (x, y) „ 
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3. Digital Filtering , A digital computer can, be used as an 
efficient tool in picture processing by using the well known techniques 
of digital filtering [2]. In order to do this, it'is first necessary to 
obtain image data in digital form. This is done by scanning the picture 
with a scanning microdensitometer, sampling and digitizing the densito- 
meter records, and thereby producing a sequence |s^ * (i, j) [ i=l,..., 

M and j=l, ..., N where 

. 


S x * (i, j) = 


h ( v yj 


) 


(1.3) 


and (x. , y.) are the points at which the picture function S (x, y) is 
measured b^ the scanner. A digital filter is just a transformation of 


the sequence 
The sequence 
on a CRT disp 


Si * (i, j) 


which produces a new sequence-|S 2 * (i, j)}. 

S 2 * (i, j)[ can be converted back to a picture, for example, 
ay unit. 


The above steps in picture processing are summarized in Figure 1.2. 


The philosophy of software development is governed to a great 
extent by the hardware available for implementing it. Hence, a brief 
description of the facilities which represent the three main blocks of 
Figure 1.2 is given in the following subsection. 


C. Description of the Processing Facilities 

The hardware used for the digital processing applications dis- 
cussed in this report is described below. The microdensitometer shown in 
Figure 1.2 is an Optronics Photoscan P-1000. It is a rotating drum scanner 
whose apertures and scanning intervals can be set to 12.5a* , 25a* and 50 a* . 
It can be used to quantize densities between 0 and 2D or 0 and 3D in equal 
steps either as eight bit (0 to 255) or six bit (0 to 63) integers. The 
range of densities can be increased beyond 3D by using a neutral density 
filter as bias. 

The microdensitometer is currently being used on its six bit mode 
which is convenient for use on the IBM 7094 (with 7- track tapes and 36 
bit words) on which all the digital processing is implemented. The 
scanner generates six bit numbers on tape and hence each 36 bit word 
corresponds to the densities at six neighboring picture points (pixels) 
on a scan line. Since the IBM 7094 is a word oriented machine it is 
necessary to unpack each of the 36 bit words generated by the microdensito- 
meter into 6 words before performing any arithmetic operations on it. 

Also, at present, the IBM 7094 has only tape units and hence all the soft- 
ware is geared towards using work tapes (not disk files) . 

There are two choices for the picture reconstruction unit in 
Figure 1.2, viz., a DICOMED 31 storage CRT display unit and an Optronics 
Photowrite P-1500. The DICOMED 31 is a 64 grey level display unit which 
converts each ,36 bit word on packed tape into six neighboring pixels on 
a horizontal line on the screen of a storage CRT, the intensity of the 
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FIGURE 1.2 STEPS IN PICTURE PROCESSING. 
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pixel being proportional to the corresponding 6 bit number. It generates 
one horizontal line per record. The screen capacity is 1024 x 1024 pixels. 
It is required that each record on tape should have at least 171 words 
(i.e., 1026 pixels). If a record has more than 171 words all but the first 
171 words of the record will be ignored. The writing on the screen is 
stopped when 1024 records are written or an end of file mark is encounter- 
ed on the tape. The Photowrite is a rotating drum device which converts 
digital data into a picture on a film. The writing aperture and interval 
can be set to 12.5// , 25// and 50// . It can convert either 6 bit and 8 bit 
numbers into pixel density on film. It is currently being used in its 
6 bit mode to be compatible with the rest of the equipment. 

The rest of this report presents the techniques and software which 
have been developed for digital image enhancement using the above facilities 
and shows examples illustrating the application of these techniques to 
radiographic images. 


D. Report Outline 

This report is presented in two main parts. The first part is 
mainly a description of the techniques and the second part is a documenta- 
tion of the software developed. 

The techniques and software are classified for convenience as Image 
Handling and Image Processing. The term Image Handling is used to cover 
such operations as unpacking and packing of the data, geometrical operations 
on image data such as transposition of arrays or rotations and translations. 
The term Image Processing refers to modification of image data arrays 
using point or local operations. Section II considers Image Handling 
operations for which software was developed, as and when the particular 
operations were found necessary. Section III describes the image process- 
ing operations, considering in particular techniques for the design and 
implementation of digital filters both in the spatial frequency domain 
and the spatial domain. Section IV presents the results of the process- 
ing operations applied to sample radiographs. 
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SECTION II. IMAGE HANDLING OPERATIONS 


In the course of developing techniques end software for radiographic 
image processing, it was found necessary to produce software for various 
image handling operations. The software, while designed mainly for radio- 
graphic image processing, includes in some cases, features which are use- 
ful for other purposes also where the inclusion of these features did not 
involve much additional labor. 


A. Automatic Picture Extraction 

The first step one has to go through in processing picture data 
after a tape of digital data is generated by the microdensitometer is to 
convert the data into a form suitable for handling on the computer. Also, 
in the interest of minimizing processing time it is desirable to extract 
as small a rectangular array of the digitized data as possible containing 
the region of interest in the picture. Further, as noted in section I-C 
it is necessary to unpack each word of digitized data into six words 
before 'performing arithmetic operations on IBM 7094. 

.. Manually, one can display the data on the tape from the micro - 
densitometer on the storage CRT screen and superimpose a grid on it and 
then determine the coordinate bounds of the region of interest relative 
to the top left corner of the screen. This, is a useful procedure when the 
picture has sufficient contrast so that?sonie features in the region of 
interest can be easily identified from the display on tjie screen. However, 
in the case of most of the structural radiographs, it has been found that 
when the original digitized version is displayed on the screen the contrast 
is so low that it is difficult to identify the region of interest. Hence 
the following technique has been developed to determine' the region of 
interest and extract it in unpacked form onto a, tape. The technique will 
be first described and the assumptions needed and its limitations are 
discussed later. •" '• . . . - ■ - v ■„>. 

In this technique, the first step is to mark the area of interest 
by masking the remaining areas of the radiograph with an opaque material 
so that the area of interest is an approximately rectangular window. Then, 
the radiograph is digitized so that the masked regions of it are included 
in the scan on all sides of the region of interest. Next, each record of 
the digitized data is read and examined for points with density numbers 
less than 63 (i.e. ,' the highest density number with, 6 bit digitization). 

The first record which has density numbers less than 63 is noted. Start- 
ing with this record the points A^ and B^ are determined for the ith 
record as shown in Figure 2.1. The point A^ is the first point in the 
ith record which has density number less than 63 and the point B^ is the 
last point in the ith record whose density number is less than 63. Also, 
the last record of interest is determine by now finding the first record 
all of whose density numbers are equal to 63. In the ideal case where 
the mask has an exactly rectangular window whose sides are parallel and 
perpendicular to the scanning direction, and where all the density numbers 
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in the region of interest are less than 63, this method will yield the 
exact boundary of the region of interest and the pixel numbers ©f A-^ 
and Bi in the ith record will be constant with respect to i. If however, it 
is known that most of the density numbers in the region of interest 
are less than 63 and that there are no long sequences of consecutive 
points with density numbers equal to 63 along a scan line at either end 
of the region (i.e., no highly dense straight lines in the scanning 
direction) then this method will yield the locations of Ai and in 
the various records very close to, but not exactly on, the beginning and 
end of the boundary of the window in the mask. Further, margins can be 
allowed on all four sides to compensate for the window in the mask not 
being exactly a rectangle and for possible skewness in its orientation 
with respect to the scanning direction. Then, to ensure that the area 
which is finally extracted is within the window marked by the mask, the 
right most location of and the left most location of B^ are chosen 
after allowing the specified vertical and horizontal margins. Thus the 
area to be extracted is determined. After finding the area to be extract- 
ed the data tape is rewound and the required part of the data is unpacked 
and written on another tape. 

This method works under the following assumptions: 

(a) The region of interest can be isolated using an 
approximately rectangular window in a mask. 

(b) The skewness in the sides of the mask is limited 
such that the right most location of A.^ is to the 
left of the left most location of B^ (after vertical 
and horizontal margins are allowed) . 

(c) The region of interest does not contain long dense 
strips in the scanning direction which give density 
numbers equal to 63 at either end of the region. 

It has been found that in most of the radiographic work done so far, 
the above assumptions are satisfied. A margin of 20 to 30 pixels on all 
sides has been used successfully in extracting the region of interest in 
several radiographs digitized with a 12. 5 /J scanning interval. 

The details of the program APES (Automatic Picture Extraction and 
Scaling) which implements the above technique are presented in Part II 
of this report. This program includes, as options, linear stretching of 
density numbers between 0 and 63 on the basis of the maximum and minimum 
density numbers in the region of interest and overriding the automatic 
determination of the region to be extracted and extracting an externally 
specified rectangular region. 


B. Transposition and 90° Rotations 

In many picture processing applications, it is necessary to perform 
filtering operations in two dimensions. Several’ two. dimensional operations 
on pictures lend themselves to implementation as two one dimensional opera- 
tions either in tandem or in parallel. Such operations can be referred to 
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as separable operations and are much simpler to implement than general two 
dimensional operations which are not separable . 

One dimensional operations handle one record of data at a time,, 
Therefore, a horizontal filtering operation (i.e., along the scanning 
direction) on digitized picture data can be performed directly from the 
tape generated by the picture extraction software described in section 
II-A. However, to perform vertical filtering (i.e., perpendicular to 
the scanning direction) it is necessary to have one point from each of 
the records on tape. Therefore to simplify data handling it is necessary 
to transpose (or rotate through 90°) the data array so that pixels along 
a straight line perpendicular to the scanning direction will form a record 
on tape. Since the data arrays are almost always too large to fit in core, 
it is necessary to use a special technique to transpose (or rotate) the 
data arrays. This technique is described below. 

This technique treats the given matrix A as made up of several 
128 x 128 submatrices. If the number of rows or columns of A is not a 
multiple of 128, the matrix A is padded at the ends with arbitrary 
numbers which are immaterial so that it can be treated as a composite of 
128 x 128 submatrices. Denote the 128 x 128 submatrices of A (padded if 
necessary) by A^, A i2» •••> A ln> •••> •••> Ami’ •••» A^, as shown in 
Figure 2.2. Also, let Aj_, A 2 , ..., A n denote the submatrices of A consist- 
ing of the first, second, ..., nth 128 columns of A. 

Step 1: The program reads the input matrix A row by row and splits 

each row into n parts and writes each part on a work tape. Thus, at the 
end of one reading of the input tape, the ith work tape will have the sub- 
matrix Aj_ for i = 1, ..., n where each row of the matrix is written as 
one record. The input tape and work tapes are rewound. 

Step 2: The submatrix An is read from the first work tape (i.e., 

the first 128 records on the tape) and A^ is written as one logical 
record on the output tape. A 12 is read from the second work tape and AJ 2 
is written as the second logical record on the output tape. Thus, at the 
end of this step the output tape will contain mn records, in the order 
A ll’ A 12’ •••’ A ?n’ •••> •••» A-£i, ..., A^. The output tape and the work 
tapes are rewound. 

Step 3: Now, Aj^ is read from the output tape as one logical record 

and written on the first work tape as 128 logical records. Similarly, 

A^ 2 > • ••> A? are written on the first work tape. Thus, aT^, . .., At are 
written on tne i fc h work tape for i = 1, ..., m. The work tapes and tne 
output tape are rewound. 

Step 4: Now, note that the first logical record on the i*^ work 

tape is the first row of Ajj. Therefore, the first record is read from 

each of the work tapes, and the first row of A^ is formed and written on 

the output tape as one logical record, remembering to throw away the 
arbitrary numbers, if any, in A^ which were used for padding. Thus, all 
the N rows of AT are written on the output tape as one logical record each. 


11 



A = 



-An 

*» 


1 A t “I 

In 


A 21 

A 22 

I — — 


r | 1 | I 

— 

— 

| 

| 

= [Al | A 2 | . . . | A n J 

• 

• 

• 

• 

l * 

1 : 


— A ml 

A m2 

| • • • 

1 A m _ 


_a t 

An 

a t 

1 A 21 

1 

| • • • 

j k h~ 


a t 

A 12 

1 a t 
| A 22 

l * * * 



• 

• 

• 

| • 
• 

| • 

| • 

• 

| • 



1 aT 
t 2n 

| * * * 

1 ftI 
1 rnn — 



Figure 2.2. 


The Partitions of the Matrix A and its Transpose 
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The above method needs only a very small modification in order to 
generate 90° rotations instead of the transpose of the data array. 

First of all, note that if ^ denotes the rotation of an M x N 
matrix A about its center by 90° in the counterclockwise direction, then 

a ij = a i,N-i+l for i=l,...,N 

j=l,...,M (2.1) 


where Xij denotes the i tb row, j tb column element of the matrix X. 

V 

Also, if A denotes the rotation of A about its center by 90° in the 
clockwise direction, then it follows that 


v 

a ij a M-j+l, i 


for 


Therefore, it is obvious that 

A _ T 
a ij “ a N-i+l>J 

and 


i=l,...,N 
J 1 , • . . ,M 


a i j== a 1 , M-j+1 

for i=l,...,N and j=l,...,M. Hence, to generated the first step is 
modified so that when each row of A is read, the elements are first 
reflected about the center of the row and then the record is split into 
n parts and written on work tapes. To generate X, the last step is 
modified so that each record of A^ is reflected about the center of the 
record just before writing it on the output tape. 


( 2 . 2 ) 


(2.3a) 


( 2 ^ 3b) 


C. Reflection 


In many filtering operations it is necessary to find the convolution 
of an array of filter weights with the data array. If the data array is 
an M x N matrix A and the filter array is an m x n matrix G, the output 
array B is given by 

m n 

b ij = £ g rs a i-r> j-s 

r=l s=l 


(2.4) 


It is seen that the formula (2.4) gives only a (M-m) x (N-n) array 
B if only the points a^ in the array A with subscripts i=l,...,M and j=l, 

. . . ,N are to be used. In order to get a filtered picture of the same size 
as the input picture it is necessary to define a^j outside the given range 
i=i,...,M and j=l,...,N.- .The definition of a^j outside the range determines 
the values of bij in equation (2.4) for i=l,...,m and j=l,...,n. Even 
though a^j can be defined arbitrarily outside the range of the input array 
it has been found experimentally that in most cases the edges of the out- 
put picture will appear better if a^. are defined for points outside the 
range by reflecting the given picture array in the corresponding edges. 
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That is, for example, 


a ij “ a 2 -i , j 

for 

2-M < i < 1 

(2.5a) 

and ajj = a^g.j 

for 

2-N < j < 1 

(2.5b) 


Now, horizontal filtering operations require handling the data 
sequentially record by record, each record being generally small enough 
to beheld in core. In these cases, the required reflections in vertical 
edges can be easily generated in core during computation. But in computa- 
tions where reflections of data in the horizontal edges of the picture 
are needed, as in two dimensional fast recursive filters discussed in 
Section III, it is necessary to generate data files containing the 
reflections of the data array in the first and the last records. The 
number of records to be reflected depends upon the number of rows in 
the filter array. The method for generating the reflections when the 
part of the data array to be reflected does not fit in core is outlined 
below. It is assumed that M r records are to be reflected in both the top 
and bottom edges of the array and M is the total number of input records. 

Step 1: The dimensions of a two dimensional array W are specified f 

according to the number of words per record in the given data array and 
the core capacity, the purpose being to store as many records of data as 
possible in core. Suppose the number of rows in W is m. Then the number 
of work tapes needed is the highest integer in (Mj- — 1) /m . This number 
is found. The next steps are described in relation to figure 2.3 where 
the number of work tapes is assumed to be 2 . 

Step 2: The first m records of data starting with the second record 

(constituting the subarray A-^ of the input data) are read into core and 
transferred to the first work tape. Similarly the subarray A 2 is trans- 
ferred to the second work tape. Now, the remaining (M r -2m) records con- 
stituting the subarray A 3 are read into core and written on the output 
file in reverse order (subarray A 3 ) . Next the array A 2 is read from the 
second work tape and written on the output file in reverse order (array 
A 2 ) . Next, A^ is read from tt\e first work tape and written on the output 
file in reverse order (array Ai) . 

Step 3: The input file is now rewound to the beginning of the 

file and the first (M-l^-l) input records are copied onto the output file. 

Step 4: The next m records of input constitute the subarray A 5 . 

These are read into core and copied onto both the first work tape and the 
output file. Similarly, the subarray Ag consisting of m records is copied 
onto both the second worktape and the output file. Now, the next (M r -2m) 
records are read into core and the array A 7 is formed. The last record 
of input is now copied to the output file and the arrays^ A 7 , A^ and A^ 
are written in reverse order on the output file as A 7 , Ag and A^ respectively. 


D. Packing and Centering for Image Reconstruction 

After the processing operations are performed on the image data 
and the resulting numbers are converted into integers between 0 and 63 
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FIGURE 2.3 REFLECTIONS OF DATA ARRAY A. 
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(by some type of scaling if necessary) , the data must be packed so that 
it is in a form required for the image reconstruction devices. Also, it 
is desirable to arrange the data so that the imhge when reconstructed is 
properly centered on the screen or film. Further, in many cases where it 
is necessary to study the effects of different filtering operations on 
a picture it is useful to display several images on the screen at a time. 
Therefore, a program has been developed which will prepare data for multiple 
frame display (or multiple frame writing on film). This program can also 
enlarge a picture by specified integral factors in vertical and horizontal 
directions by repetition of pixels. A brief description of its operation 
is given below. 

It is assumed that the number of frames to be packed is a composite 
number N and its factors, viz., the number of frames in the vertical (N ^ 
and horizontal (N^) directions are specified. Also, the enlargement factors 
are assumed to be given. In the case of preparing data for display on the 
screen the screen size (in the present set up) is known to be 1024 x 1024 
pixels. In the case of preparing data for writing pictures on film it is 
assumed that the film width and interval desired between vertical frames 
is specified in millimeters and the writing interval (i.e., interval bet- 
ween pixels) is specified in microns. The program first computes the 
horizontal and vertical margins between frames in terms of number of pixels. 
Next, if N^>1, the x N v frames are packed and written on work tapes, 
the ith work tape containing the i*-* 1 set of N v consecutive pictures, 
repeating pixels if enlargement is required and the areas corresponding 
to margins being filled with zeros. The i fc ^ record to be written on the 
output tape is formed by combining the ith records on each of the work tapes, 
for i=l,2,3,... . 


E. Other Image Handling Operations 

There are several other image handling operations which have been 
found useful. These operations are very simple to perform and do not warrant 
a formal description of techniques. Some of these operations and their 
purposes are discussed below. 

1. Translation of Pictures . The term translation here refers to 
shifting of picture data by specified amounts in the horizontal and vertical 
directions. This is essentially a re-indexing operation and is useful for 
registration of two images. It is well known that point by point averaging 
of several pictures of an object reduces the random noise content of the 
image and hence improves the signal to noise ratio. But in performing the 
averaging digitally, it is necessary to align (register) the data so that 
the same point on the various pictures of a scene has the same address in 
the digitized files of these pictures. Therefore, a program has been 
developed for translating the picture data to obtain proper registration 
given the amount of shift needed. 

2. Reduction of Dimensions of a Picture . When an area of a 
picture is digitized with high resolution it is possible that all of the 
image cannot be displayed on the screen. Also, it may be desirable, some- 
times, to work on an image first with a reduced resolution so that the 
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entire image can be processed quickly and an area of interest can be identi- 
fied for further processing. Therefore, a program has been developed for 
reducing the dimensions of a given picture data file by given integral 
factors either. by skipping points or by averaging density numbers over 
rectangular blocks of points. 


1 
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SECTION III. IMAGE PROCESSING OPERATIONS 


This section presents the mathematical details of the image process- 
ing software developed. This section is divided into two main subsections. 
The first subsection considers the filter design techniques from either 
frequency or spatial domain specifications. The second subsection considers 
filter implementation techniques. 

A. Filter Design Techniques 

In picture processing it is necessary to take into account the 
nature of the object whose image is being processed, and also the nature 
of the corrupting noise which might be due to various sources such as 
fluctuations in the light source and the electrical noise in the output 
of the light sensing device in the image scanning system. A useful 
technique for describing the nature of the object and the various noises 
and modifications introduced by the imaging system is the "frequency 
domain" method which is based on the Fourier Transform. This method 
is useful when all the systems and operations under consideration can be 
assumed to be linear and spatially invariant and the noise can be assumed 
to be stationary. When such assumptions are not valid one has to design 
and implement the filtering operations indirectly in the spatial domain. 

Both the frequency and spatial domain methods of filter design are dis- 
cussed below. 

1. Frequency Domain Approach: In this method, as in classical 

Control Theory and Network Analysis, the inputs (i.e. the spatial varia- 
tions of the intensity of the light from an object) to the imaging system, 
the outputs from the imaging system (i.e. the spatial variations of the 
intensities recorded on the film) and the imaging system itself are 
described in the Fourier transform domain. In other words, the inputs 
and outputs are decomposed into sinusoidal components of varying frequency 
amplitude and phase and the imaging system is described by its "Optical 
Transfer Function" (OTF) which is the relative amplitude and phase response 
of the system to sinusoidal inputs of varying frequency. If the imaging 
system is linear, it is well known that the output spectrum (i.e. the 
Fourier transform of the output) is just the product of the OTF of the 
system and the input spectrum. [3] 

Also, if the input spectrum is band limited then the output spectrum 
is also band limited. Thus, if we are interested in finite resolution (as 
is most commonly the case) then the spectrum of the pictures produced by 
the imaging system would be band limited. Then, if the picture is digitized 
using a sampling rate greater than twice the maximum frequency present, 
the spectrum of the digitized picture is the same as that of the original 
picture over the frequency range of interest. Therefore, even the digital 
components of the image processing system such as the digital filter can 
be described using the frequency domain approach, and the output spectrum 
of any component is the product of the input spectrum and the transfer 
function of the component. 
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In order to apply the concepts of frequency domain analysis used in 
control theory to the analysis and design of image processing systems, it 
is necessary to use a two dimensional Fourier transform instead of the 
usual one dimensional Fourier transform since there are two independent 
(spatial) variables (x, y) in place of one independent (time) variable t. 

For simplicity, the analysis will be carried out using one dimensional 
Fourier transform and the analysis can be easily extended to two dimensional 
filtering. 


a. Choice of filter transfer function: The first step in design- 

ing a digital filter in frequency domain is to obtain the specifications of 
the filter. In the frequency domain, specifying a filter consists in giving 
the desired range of frequencies that are to be passed unchanged and the 
range of frequencies that are to be stopped, attenuated or amplified. For 
example, a hairline crack in a metal or a bone is a high frequency phenomenon 
compared to the rest of the object and hence if the crack is to be made 
prominent the low frequencies must be suppresssed. Thus, a high pass filter 
should be used. Some of the other classes of filters are low pass, band 
pass and high emphasis. 


The magnitude of an ideal low pass filter transfer function is shown 
in figure 3.1a. The transmission in the pass band is unity and that in the 
stop band is zero. There is a discontinuity at the transition between the 
bands. Transfer functions can be digitally realized with fewer terms (and 
hence yield faster filter implementations) if they are continuous functions 
of frequency [2], Also, sometimes, the desired transfer functions might be 
such that there is a finite gap between the pass band and stop band which 
may be called a "don't care" band where the transfer function can be 
arbitrary. Therefore, it is desirable to approximate an ideal filter by a 
continuous transfer function. A simple continuous approximation to the 
ideal low pass filter is the n c ^ order Butterworth filter defined by 

G (w) = [l + (w/« Q ) 2n j" * (3.1) 

Where n is an integer and oj q is the cut off frequency. The accuracy of 
the approximation to the ideal filter increases as n, the "order" of the 
filter increases. 


Another type of approximation that can be used when there is a "don't 
care" band is [4] : 


G (w) = / 0 
1 


for 


>«> 

— T 


(3.2) 


for ai < w 


. w 

L\ p 

for w < I <jj 

W C ) 

c — 1 




where p is a positive integer and[ U( -,, w ] is the "don't care" band. 
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FIGURE 3.1a IDEAL LOW PASS CHARACTERISTIC. 



FIGURE 3.1b BUTTERWORTH APPROXIMATION. 



FIGURE 3.1c FREQUENCY RESPONSE DEFINED BY EQUATION 3.2 



FIGURE 3dd FREQUENCY RESPONSE DEFINED BY EQUATION 3.4 
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The function G (w) given by (3.2) is continuous but does not have a 
continuous derivative. The following modification of G (w) assures 
continuity of the first derivative for p >2: 

G (<*>)= / 0 for | oj | 2 w X 

1 1 1 for \<*j\< 0J C 


(OJ T - )P 


(OJ i£_ OJq) 

i. (M-" c ) p 


for co 


o-l <_\0) 


(p-1) 


<t\ - w c ) ( w ! - w c ) *' for < 

where is an arbitrary number such that uiq <w^ <w^. 

In particular, choosing = ( uq + u T )/2 we get: 

fori oj > co m 

I J— T 

f or|o>|< o) c 


OJ 


<_co T 


< OJi 



'T ~ *"c / f° r w l <J4>|< W T 

X - | I w I - U C 


l« T -W c 


for <_ W 


< CO l 


(3.3) 


(3.4) 


The forms of approximating characteristics represented by equations 
(3.1), (3.2) and (3.4) are shown in figures 3.1(b), (c)- and (d) respectively. 

Approximation to a high pass filter transfer function can be obtain- 
ed by (1 - G(w)) where G (w) is an approximation to a low pass filter transfer 
function. Approximation to a band pass filter transfer function can be obtain- 
ed as a difference between approximations to two low pass filter transfer 
functions . 


Sometimes it might be desired to compensate for the errors induced 
by effects of the imaging system. If the transfer function is known, or 
may be estimated or measured and is represented by, say, G (u>) , then, 
the required filter transfer function is H (w) = 1/G (u>) . For example, 
compensation for finite size of the scanning beam can be obtained using a 
filter transfer function shown in figure 3.2. [ 5 ] 

Since various picture processing applications need several types of 
filter transfer functions it is desirable to have a general software pack- 
age which can approximate any given filter transfer function. Such a 
package has been developed and is described in the following sections. 
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First of all, the filter input/output relations are described for 
the continuous and discrete cases. 

b. Filter equations. Consider first the equations for the 
continuous case. Let the intensity function input to the filter be given 
by S(x). Let the filter transfer function be given by G (w) and let g (x) 
be the corresponding impulse response. Then, the output S Q (x) from the 
filter is given by 

S o (x) _ 

Now, if S(w) and S c (w) are the Fourier transforms of the input and the out- 
put of the filter respectively, then it follows that 

SqCcj) = G(tu) T(tu) (3.6) 

The corresponding expressions for the discrete case are as follows: 

S on s E g k S n _ k for -co<n<co (3.7) 

k= -oo 

where |S n |, | g n | and | S Qn | are samples of the input, filter impulse response 
and the output respectively at intervals of 4x. 

If the input and the filter impulse response have spectra which are 
zero outside the interval -1/(2 /ix) <. f < l/(2 4x), then the expressions 

for the Fourier transforms have the form 

S((J) = S n exp(-jWn4x) 

n= - 0 ° 

G(W) = £ Sn exp(-jWnZlx) 
n= -oo 

in the interval - n/ A x <_U)< 7rJ A x and S(W) = G(0>) 
output spectrum is also band limited and is given by 

1T 0 (6J) = G(GJ)S(aJ) (3.10) 

The numerical errors associated with assuming a function to be band limited 
and deriving equations (3.8), (3.9) and (3.10) are discussed by Coo, ley 
et al. [6] 

The process of designing the digital filter consists in choosing the 
'weights' g n such that the desired band limited transfer function G(w) is 
obtained. However, since the picture data to be filtered is. a finite set 
of numbers and obviously an infinite summation cannot be computed on a 
computer, it is desirable to design filters with a finite set of weights 
g k , k = 0, . .., L-l. Using this set of weights, the equation (3.7) can 
be written as 


(3.8) 

/(3.9) 

/ 

/ 

/ 

= 0 elsewhere. The 


-yi 


g (u) S (x-u) du 


(3.5) 
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L-l ... 

S on = 2 Sk s n -k n = 0, N-l (3.11) 

k=0 

noting that the input and output are also finite sets of numbers. 

v \ 

Now, in order to obtain G(w) exactly, in general, an infinite 
sequence of weights g^ is necessary. Hence, a finite sequence of weights 
will only provide an approximation to G( w) . Filters with smaller number 
of weights generally result in faster implementation but higher errors 
in approximation. There are several techniques in the literature for 
obtaining filter weights to approximate given transfer functions. Two 
of these methods will be described below. The first of these can be used 
to approximate an arbitrary continuous transfer function G(w) . Note 
that discontinuous transfer functions cannot be approximated to arbitrarily 
small errors even with an infinite sequence of filter weights [ 2 ]. 

c. Helm's 4-T's method [ 7 ] (transform-truncate-transform-test). . 

Let G(w) be the specified filter transfer function which is zero outside, 
the interval 

- 7T/ A x < 10 < tt / A x. 

(a) First, a large number N of samples is chosen such that all 
the significant features of the frequency response such as the peaks, zeros 
and discontinuities are sufficiently accurately represented by the sampled 
response and a sequence |G n | is formed such that 

G n = G(CJ n ) (3.12) 

where 

(J n = 2 TT . n for n = 0, 1, ...,(N/2)-l 

4x N 

and (J n = -2 7T N-n for n = N/2, ...» N-l (3.13) 

A x N 

(b) Now, the IDFT (Inverse Discrete Fourier Transform) of 
the sequence |G n | (n = 0,1,..., N-l) is computed. The IDFT, denoted by 

jg' n j gives the impulse response corresponding to the frequency response |G n }. 

(c) A number L<N is chosen. In the sequence •|g' n [, N-L 
contiguous values are set equal to zero, g 0 ' and g'^_^ being considered 
contiguous. The sequence |g n | is obtained by permuting ^ ' n | such that the values 
which were set to zero appear at the end. That is, g n = 0 for n=L, ...,N-1. 

(d) The DFT of the sequence |g n |is computed. Let this DFT 

be denoted by|G' n |. The error e in the approximation is defined as a suit- 
able function o f | G n - G' n |. For example 

e = Max G„ - G' (3.14) 

0 < n < N-l n 
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(e) If the error e is not small enough the steps (c) and (d) 
are repeated with a larger value of L. If e is smaller than the specified 
error then L is decreased and the steps (c) and (d) are repeated. The 
new values of L can be obtained either by a successive bisection process 
or by a modified Newton -Raphs on Method. 

This method is illustrated in figure 3.3 

d. Programs to implement the 4-T's method: A set of programs 

has been developed to implement the method described in the steps (a) through 
(e) in the previous subsection. The programs have two options. The first 
option is to double the value of L whenever the error e is greater than or 
equal to the s£ecified error. When the_ value of L exceeds a prespecfied 
maximum value L , L is taken equal to L and the impulse response is truncated 
to L terns and written as an output tape. The second option is to use 
successive bisection of intervals of L. First the error is found with a 
small initial value L-^ of L (which is less than L ). Next the error is 
formed for L = L. Now for the i fc ^ iteration the value of is given by 

Li = (Li_l + Lj) /2 for i>3 (3.15) 

where 

j = Max { j j j < i-1 and (e^ - E) (e^ - E) <^0 } (3.16) 

and 

e.^ is the error during the i th iteration and E is the given error. 

Several types of filters have been designed using these programs. Figures 
3.4, 3.5 and 3.6 show the magnitude versus frequency plots of some examples 
of the original transfer functions for low, high and band pass filters 
respectively and the approximating transfer functions produced by the above 
method. The magnitude characteristics for these filters were assumed to be 
derived by a basic low pass characteristic given by equation (3.2). The 
maximum permissible error e in the equation (3.14) was assumed to be 0.05 
(i.e. maximum gain of approximately -26 db in the stop band and an error of 
0.22 db in the pass band). The number N of samples of the frequency response 
was taken to be 256. The parameters defining the filters shown in figures 
3.4, 3.5, and 3.6 and the values of error for the attempted values of L 
are shown in tables 3.1, 3.2 and 3.3 respectively. The first option (doubling 
L until error is below the specified value) was used in these cases. 
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FIGURE 3.3 STEPS IN FILTER DESIGN (4-T METHOD) 
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FIGURE 3.3 (CONTINUED) STEPS IN FILTER DESIGN (4-T METHOD) 




TABLE 3.1 APPROXIMATIONS TO A LOW PASS FILTER 


Pass Band = [0, 0.15w s ], Stop Band = [0.25 , 0.50w s ] (The pass band and stop 
band are joined by a third order polynomial as in equation (3.2)) 


Filter Length 

Error 

L 

e 

8 

0.3265 

16 

0.1633 

32 

0.0728 

64 

0,0285 


TABLE 3.2 APPROXIMATIONS TO A HIGH PASS FILTER 

Pass Band = [o.25w s , 0.50w s ], Stop Band = [ 0, 0.15w s ] (The pass band and stop 
band are joined by a third order polynomial as in equation (3.2)) 


Filter Length 

Error 

L 

e 

8 

0.3265 

16 

0.1633 

32 

0.0728 

64 

0.0285 : 


TABLE 3.3 APPROXIMATIONS TO A BAND PASS FILTER 

Pass Band = [0.10w s , 0.25w s ], Stop Band = [ 0, 0.50w s J U [ 0.30o> s , 0.5w s ] 

(The first and second stop bands are connected to the pass band by a 2°d and 
3 rd order polynomial respectively. 


Filter Length 
L 

Error 

e 

8 

0.3766. 

16 

0.3000 

32 

0.1820 , 

64 

0.0895 

128 

0.0371 
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Figure 3.4 Low Pass Filter 

(a).' Approximating Frequency Response 


(b) Specified Frequency Response 
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Figure 3.5 High Pass Filter 

(a) Approximating Frequency Response 

(b) Specified Frequency Response 
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e. Frequency sampling method of design Helm's 4-T's method 

described above first finds the filter impulse response and truncates it 
to a desired degree of approximation in order to obtain a short filter. 

Another technique for approximating a class of filter transfer functions is 
due to Rabiner et al [8] who uses a frequency sampling approach. In this 
method the resulting impulse response is ’folded' or 'aliased' rather than 
being truncated. The steps in the design procedure are summarized below. 

(a) A set of frequencies is chosen at which the sampled 
frequency response is specified. This set consists of a number of points 
much smaller than that required for an accurate representation of the 
desired frequency response. The values of the transfer function at some 
of the frequencies in the above set are left as design parameters. For 
example, in the design of an ideal low pass filter the transfer function 
at the pass band frequencies is fixed at 1 and at the stop band 
frequencies is fixed at 0, but a few frequencies are considered transition 
frequencies joining the pass and stop bands and the values of the transfer 
function at these frequencies are considered design parameters. 

• . .’‘ft 

(b) The values of the continuous frequency response are obtained 
in terms of assumed values of the design parameters by using either an explicit 
formula (given by a sampling theorem) or interpolation using fast Fourier 
transform or chirp z-transform algorithm. Here, continuous frequency response 
means the response at a number of frequencies which is much larger than the 
number of frequency samples chosen in step (a) and adequate to characterize 
the necessary features of the desired transfer function. 

(c) * Using the interpolated frequency response the values of the 
design parameters are readjusted until an error criterion is satisfied. 

The error criterion, for example, could be that the maximum of the absolute 
value of the difference between the approximating and the desired frequency 
responses in a given frequency range be a minimum or in the case of an ideal 
low pass filter that the maximum gain in the stop band should be a minimum. 

(d) When the error criterion has been satisfied, the final values 
of the free parameters are used along with the fixed values of sampled transfer 
function which were determined in step (a). 

Rabiner et al have used this method to design a large number 
of low pass and band pass filters and wide band differentiators. They have 
provided tables in reference^] giving design data. These tables give the 
band width, total number of frequency samples used, the values of transition 
samples and the minimum gain in the stop band achieved for the low pass and 
band pass filters. High pass filter frequency samples can be obtained from 
these tables by subtracting all the low pass filter samples from 1. The 
frequency samples for band widths not listed in the tables can be obtained 
by linear interpolation of the tabulated data. It has been found experi- 
mentally [ 8 ] that the deviation of the maximum gain in the stop band obtained 
by linear interpolation will be less than 6 db from the optimum. 
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2 . Spatial Domain Approach : The main difference between the 

frequency domain approach to filter design and the spatial domain approach 
occurs in the way the desired characteristics of the filter are specified. 

In cases where the signal and noise spectra are known and are approximately 
non-overlapping it is simple to choose a filter transfer function. This 
would be a frequency domain approach. On the other hand, in some cases such 
as matched filter design it is much simpler to specify the filters directly 
in terms of the weights in the spatial domain. Also, nonlinear and spatially 
variant filtering operations cannot be specified conveniently in the frequency 
domain. This subsection presents two types of filters which are designed 
using the spatial domain approach. 

a. Matched filters (template matching): Matched filters have 

traditionally been used for detection of signals in presence of noise in 
communication systems [9-11]. These are derived in communication theory for 
the one dimensional case and are optimal in the sense that they maximize 
the signal to noise ratio. A derivation of the matched filter for the two 
dimensional case which is useful for picture processing is given below 
for the sake of completeness. (See reference [l2] for a more detailed 
discussion). The continuous case is considered first and the c orres ponding 
equations for the discrete case follow easily. 


Let fj (x,y) be the input signal and let the corrupting 
noise be n (x,y) . Assume that n (x,y) is additive noise generated by a 
wide sense stationary random process which is ergodic in its autocorrela- 
tion. Let Rjj (T,T) be the auto correlation and S^i (u,v) be the power 
spectrum of noise. Let h (x,y) be the impulse response of the filter. 


Now, the input to the filter is given by 

f (x,y) = fj. (x,y) 4- n (x,y) (3.17) 


The output of the filter due to the signal is given by 
f Q (x,y) = fj (x,y) * h(x,y) where * denotes convolution. Therefore, the 
output signal has instantaneous energy at the point (£ ,TJ) given by 

s -1 £ 0 i I 2 




(u, v) exp | j (u | + V*7) ] du dv 


I If Fj (u,v) H (u,v) exp j j (u£+v7?) } du dv 


(3.18) 


where Fj, F Q and H are the Fourier transforms of fj, f Q and h respectively. 
The power spectrum of the output of the filter due to noise input is given by 


^NO ( u > v ) ~ ( u > v ) 


H (u,v) 


2 


(3.19) 
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Therefore, the signal to noise ratio of the output is given by 


J 

f A (u,v) H (u,v) exp {j (u$ +V77)} du dv 

2 

1 J 

'J 

fs ni 0*,v) 

H (u,v) | 2 du dvj 



(3,20) 


Now, it' is well known (see [3] , for example) that the power spectrum S™.(u,v) 
is non-negative for all u,v. In particular, assume that Sjqi(u,v) > 0 for all 
u,v. (This is certainly true for the case of white noise input. For a 
general treatment of the case where S NI (u,v) is permitted to be zero for some 
values of u,v see [13] ). Then, using the Cauchy-Schwarz inequality on the 
numerator of the right side of equation (3.20), it follows that 


S < 

N “ 


/ /|Fj (u,v) 


/S NI (u , v) dudv 


f Ac- 


V) 


H(u, v) 


^ dudv 


J J S NI < u ’ 


v) H(u, v) 


dudv 


(3.21) 


That is. 


| - If | F i 


v) 


/S N i (u,v) du dv 


(3.22) 


From (3.20) and (3.22) it is easy to show that the maximum value of the output 
signal to noise ratio is achieved if and only if 

H (u,v) =XFj (u,v) exp (-j (u£+ V 77 ))/S NI (u,v) (3.23) 

which is the condition for (3.22) to hold with equality, where X is a constant 
with respect to u and v and the asterisk denotes complex conjugation. If the 
noise is assumed to be white, then the expression (3.23) reduces to 

H (u, v) = Af* (u,v) exp (-j (u£+v77)) • (3.24) 

and the impulse response (i.e, point spread function) of the filter can be 
written as 

h(x,y) = Xf ! ( £-x, 77-y) (3.25) 

Therefore, for the particular case of additive white noise in the input, the 
output signal to noise ratio at the point (£,7)) is maximized if the filter 
point spread function '’matches" the input signal shifted to the point (£,*7) 
and reflected point by point in (£,?7) to within a multiplicative constant. 

In other words, if the location or locations of a particular image pattern 
is to be detected in the presence of additive white noise one would just move 
a template having that pattern over all points of the given picture and find 
the locations at which the maxima of the cross correlation between the given 
image and the template occur. Further, if contrast variations within a picture 
need to be compensated for, then the maxima of normalized cross correlation 
should be examined. 
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The discrete version of the matched filter equations can be 
written down immediately. We shall assume here that only the variations 
of the picture function over the local average are significant, so that 
the template can be assumed to have zero mean. A set of weights ^g^ 
are chosen with (k, £)€[-K, K] X [-L ,L] so that they represent the shape of 
the pattern to be detected. The input picture (i.e., the picture to be 
filtered) is given hylx^j}® {s^j + n.jj} where n-y is white noise which 
is ergodic in its mean and variance ana statistically independent of Sjj. 
Also, n^j has zero mean. The normalized cross correlation between {gf 
and | x | is given by 


yij = E 8 ke *i+k,j+£/( 
k,£ 


iji 


) 


(3.26) 


where 


II i 


and 


Is II- (£« l O' 

Ihijil = (E 

k,£ 

Now, the expected value of yij is maximum when 8i+k, j+ t~ Ag^ for 
K== -K,...,K and £= -L,...,L. Further, due to the assumptions on the 
noise, the actual values of y^j will be approximately equal to the expected 
value for sufficiently large values of K and L and is given by 


(3.27a) 

(3.27b) 


¥ 

yij ~ 


J°kl i+k. Vrl 

|;«|| <1 s ijf ^ + | n ij 




(3.28) 


When a 'match' occurs, -the value of yij is given by 


'ij « 


<l|8i| Z + 


n 




(3.29) 


noting that || n i j if is approximately constant with respect to i and j and 
denoting it by j|n|| . Thus, if one knows the statistics of noise, one can 
use equation (3.29) to obtain a threshold on the normalized cross correla- 
tion to decide whether the template matches at a particular location. 


In recent years, the application of matched filtering for optical 
character recognition and image pattern recognition has become very common 
[l, 12, 14, 15 }. Unlike other image enhancement techniques matched filter- 
ing does not, in general, produce an output image which looks like the 
pattern being enhanced. Instead, it produces a "correlogram" from which 
only the locations of template match can be determined. One can then 
reconstruct the desired image using the correlogram and the template. 

In some cases, however, the correlogram itself (appropriately thresholded) 
is a reasonable representation of the desired image. A particular case 
where this is true is when lines of a specified width in a particub r 
direction are to be enhanced. t 
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b. Density Manipulations: Even though all filtering operations 

described so far are in fact manipulations of the densities in the input 
picture to get an output picture, we shall specialize the term density 
manipulations to mean point operations (in contrast with local operations) 
where the output due to a particular picture point depends on the density 
at the point only (and not on the densities at some of the neighboring 
points also). Such operations are sometimes useful for image enhancement 
and are much faster than local operations. 

(1) Table look up method: When the input data consists of only 

a small set of integer values compared to the total number of data points 
the density transformation can be implemented by defining it as a table. 

When the input data is read the output is found by merely looking up the 
table. This is a very fast method since no arithmetic operations are 
involved. Examples of this method are (i) contrast stretching by subtract- 
ing the minimum density number from the density number at a point and 
linearly rescaling the density numbers over the available range (0 to 63) 
(ii) modifying the distribution of densities by remapping of densities, say, 
to yield a linear distribution curve. 


(2) Linear and non-linear scaling: More often than not, when image 

data is filtered, the resulting set of numbers are both positive and negative 
and are in floating point format (particularly in high pass and band pass 
filtering). Therefore, in order to be able to display the filtered image 
it is necessary to convert the numbers to integers between 0 and 63. In 
order to obtain as much contrast as possible in the displayed image we set 
the smallest number in the array to be 0 and the largest to be 63 and rescale 
the data using a mono tonic and non -decreasing function. The mono tonic 
function could be either linear or nonlinear (square law, logarithmic, 
exponential, etc.) depending on whether contrast stretch should be uniform 
over the density range or whether contrasts in some density ranges are to 
be emphasised to a greater extent than in others. Sometimes it might be 
desirable to set all data below a certain value to 63 and scale the data 
in between according to a monotonic scaling function. 

A program "SCALE" has been designed to take all these possibilities 
into account. This program handles the input-output operations, computes 
output using an externally specified scaling function, rounds off the computed 
values to integers, truncates the resulting set of integers (if necessary) 
to lie between 0 and 63 and, as an option, computes, the maximum, minimum and 
mean of the input file which can be used in the scaling function. This 
program can handle several input files on tape at a time without having to 
rewind the tape several times. (The input tapes will be rewound only once 
if the option to compute maxima, minima and means is used). 


An example of a function which can be used as an external function 
in the routine SCALE is given below: 


y ij = «i + a 2 { < x ij-°9 ^ / < of 4- o y f p 


(3.30) 
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where.o*^ is a bias parameter, a 2 >0 is the contrast factor, 0*3 is the back- 
ground to be subtracted and Ct 4 can be taken, as 


0*4 = Max ( 


Max x 




Min xij 


) 


(3.31) 


the maximum and minimum being over ij taking on all values within the domain 
of data. 


B. Filter Implementation Techniques 

A general approach to linear shift invariant filter implementation 
would be to store the two dimensional filter (either the impulse response 
or the frequency response) in core and find the convolution of the data with 
the filter by bringing into core the appropriate sections of data. This is 
the approach taken in most installations where special purpose hardware is 
available for computing convolutions and disk files are available for inter- 
mediate storage and retrieval of data. However, on a tape based system it 
is simpler to do filtering in one dimension so that the data from tape may 
be read, filtered and written on output tape sequentially, one record (i.e. 
one picture line) at a time. Two dimensional filtering may then be per- 
formed as a two pass procedure. This restricts the class of filters which 
can be implemented to 'separable' filters whose two dimensional transfer 
function H(u,v) can be written as a product Hi(u)H2(v) or as a sum 
Hi(u) + H£(v) o The two pass procedures for the product and sum type separable 
filters are shown in figures 3.7 and 3.8 respectively. 

There are several methods of implementation of one dimensional digital 
filters. These can be broadly classified as nonrecursive and recursive 
implementations. Both the implementations are discussed below. 

1. Nonrecursive Implementations : The output of a one dimensional 

linear filter is expressed in nonrecursive realizations by the formula 

L-l 

Yi ■ £ g k x i-k (3.32) 

k=0 k 

The values of y^ in (3.32) can be computed directly by multiplying 
the terms of { g } and | x } term by term and adding them. This would involve L 
multiplications and (L-l) additions for every point of output. The number 
of arithmetic operations can be reduced considerably in the case of a filter 
whose length is much shorter than the length of an input data record by using 
the Fast Fourier Transform (FFT) algorithm. A method which uses the FFT is 
due to Stockham £l6]and Helms [l^Jand is called the "select-saving" method. 

The steps in this method are summarized below: 

(a) A number L' is chosen where L'>L. The choice of an optimum 
L' is governed by the radix of the FFT algorithm which determines the total 
number of arithmetic operations required. 

(b) The sequence j g | is appended with L'-L zeros at the end. 

Denote the new sequence by {g'f . Also let the first L' terms of the sequence 
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FIGURE 3.7 TWO PASS PROCEDURE FOR PRODUCT TYPE FILTERS 
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FIGURE 3.8 TWO PASS PROCEDURE FOR SUM TYPE FILTERS 
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| x } be denoted by { x }. 

(c) The Discrete Fourier Transforms (DFT's) of |g'| and] x | 
are computed using FFT. 

(d) The IDFT { z f of the product of the DFT's above is computed. 

(e) It can be shown that due to the cyclic nature of the DFT 

only the values z ffl for m = L-l, L,..., L'-l are the valid values of convolution. 
Therefore only those values are retained. 

(f) Now, { x \ is redefined by discarding its first L'-Lfl values, 
reducing the subscripts of the remaining L-l values by L'-Lfl and appending 
the next L'-Lf-l values of the sequence |x } to be original L-l values ofjx } , 

(g) The steps (a) through (f) are repeated until all the values 
of | x | have been used. 

The number of arithmetic operations required for computing the 
convolution for N points is approximately 3K additions and 2K multiplications 
where 


K = 2NL'(log 2 L')/(L'-L+1) (3.33) 

This shows for long filters it is advantageous to use the select-saving method 
and the advantage over direct convolutions increases as the length of the 
filter increases. Now, the number of output values obtained with the select- 
saving method is equal to (N-L+l) when filtering a record of length N with 
a filter with L weights. The output obtained is given by the equations 

y 0 = x O g L-l + x l g L-2 + ••• + X L-1 §0 


y N-L = x N-L 8 L-1 + X N-L+l g L-2 + *** + *N-l g 0 


(3.34) 


However, it is desirable to obtain an output array consisting of the 
same number of values as the input array so that the filtered picture and 
the original picture have the same size. To do this, it is necessary to 
append the input array suitably with (L-l) extra components (which may 
be chosen either zero or some other numbers, essentially to make the filtered 
picture look reasonable even at the edges) . Also it is convenient to store 
the filter weights generated by the 4-T's method (Section III A.l.c) approxi- 
mating the desired transfer function of the filter to within a known linear 
phase difference, i.e., to within a known shift in the spatial domain. In 
order to take these factors into account, two new arrays x^ and y^ are 
defined as shown below and select-saving method is used with x' and g as 
input and y' as output. 
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Define the sequence h = 11 q, . .., hj ^-1 as the set of weights 

giving the desired approximations to the filter transfer function found 
by Helm's 4-T's method, where h i = ... = = 0. Then, these 

weights are stored by writing i, L and the L "consecutive" non-zero 
weights { g } given by 

§0 = h i+N x -L 

§1 “ ^i+Ni-X+1 


g L-i-l h Ni“l 
SL-i = h 0 
s L-i+l = h l 


(3.35) 


g L-l = h i-l 

The output will be free from shift if 

y 0 = ...+h i x_ 1 + h 0 XQ + h Nl _ 1 x : + ... 


y^ ...+h^xj c _^ -f IiqX^ + +... 


(3.36) 


y N -l = ••• + Vn_2 + Vn -1 + h Nl-l X M + "• 

where x^ are to be suitably defined for k ^ [0, N-l ] „ 

Now, the filter output due to input | x^ j using select-saving method 
is given by 

?k " *£ §L-1 + *lUl g L-2 + *•* + x kfL-lg 0 (3 “ 37) 
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for 


k = 0, 1, . . . , N-l. 


Comparing (3.35), (3.36) and (3.37), it follows that the output in 
(3.36) is obtained if 


i-2 


i-1 


-(i-1) 


To be defined 


i+N-2 


’‘N-l 


(3.38) 


i+N-1 


N 


To be defined 


N+L-2 


(NbL-2) -(i-1) 


and 


y k = y£ f °* k= 0, 1, ...» N-l 


The program which implements the above nonrecursive filtering 
technique provides for three types of definition for when k f Co , N-l ] 


(i) periodic repitition with period N, i.e., x_i = x n-1» 
x^j = Xq, etc. 

(ii) reflection about end points, i.e., x_j = x^, Xjj = Xjj_ 2> etc » 

(iii) zero outside the given range, i.e., x_^ =0, x^ ® 0, etc. 

This filter.'implementation was tested by filtering a square wave record 
consisting of. 256 points using a high pass filter. The input record consisted 
of four consecutive, 5 values of 22 followed by four values of 18 and so on. 

The high pass filter; used is shown in figure 3.5. The original data and the 
filtered data are shown in figures 3.9, 3.10, and 3.11 respectively for the 
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three methods of appending the edges mentioned above. 


The DFT (256 point) of the original data shown in Figures 3.9, 3.10 
and 3.11 is given by 

X = 5120- X = 256 = x* X = = X* 

0 ’ 32 l-exp(j7r/4) 224 96 l-exp( j37r/4) 160 


and 

*n " 0 


for all other n € [ 0, 255 ] . 


Now, the values of the filter transfer function corresponding to nonzero 
values of the input spectrum are given by (see Figure 3.5). 


G 0 - ° : G 32 * °224 - ° ; °96 - <60 


Therefore, the DFT of the output is given by 


Y 96 - Y 160 


256 


1 - exp ( j3 7T /4) 


1 


and Y n = 0 for all other n € [0, 255 ] . 

Taking the IDFT of { Y n | , it is seen that 

0 1 i + JIT 2 l J 

The values of y n derived above are correct for all n€ [o, 255] under the 
assumption that the sequence | is periodic with period 256. This can be 
seen in Figure 3.9 which corresponds to making the data periodic with period 
256. 


However, in Figures 3.10, and 3.11, the values of y n derived above 
are valid only in the middle regions where the 'edge' effects are absent. 
Note that the edge effect is more pronounced in Figure 3.11 than in Figure 
3,10. This happens because introducing zeros at the ends of the data 
sequence causes large jumps at the ends while reflecting the data about the 
ends does not produce such severe jumps. 

This nonrecursive filter implementation was also used on a test 
pattern made up of density variations along a record according to square 
waves at various frequencies. The pictures of the test pattern before and 
after filtering are shown in figures 3.12, 3.13, 3.14 and 3.15. The low, 
high and band pass filters used in these pictures are the ones shown in 
figures 3.4, 3.5 and 3.6 respectively. Note that the filtering is only one 
dimensional. The low pass filter has a smoothing effect which worsens the 
contrast of the original picture and the high pass filter produces contrast 
enhancement. The band pass filter enhances the contrast, but not as 
effectively as the high pass filter, since the highest of the frequencies 
which contribute to the sharpness of the edges have been filtered out. 

Note that in the high and band pass filtered pictures long horizontal edges 
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Figure 3.9 High Pass Filtering of a Test Pattern 

(edges appended by periodic repetition) 

(a) Filtered Signal (b) Original Signal 
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Figure 3,10 High Pass Filtering of a Test Pattern 
(edges appended by reflection) 


(a) Filtered Signal 
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(b) Original Signal 
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Figure 3.12 A Test Pattern Generated with Square Waves 
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Figure 3„13 Low Pass Filtered Test Pattern 
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Figure 3.14 High Pass Filtered Test Pattern 
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Figure 3.15 Band Pass Filtered Test Pattern 
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are missing since they are of zero frequency in the direction in which 
filtering is carried out. 

Also note the absence of the vertical lines corresponding to the 
highest frequency in the low and band pass filtered pictures (the top and 
bottom l/16th s Q f the pictures). 

Figure 3.16 shows the effect of band emphasis filtering on the test 
pattern of Figure 3.12. The filter transfer function used is the same as 
in Figure 3.6 except that the stop band gain is changed to 0.12 from 0. In 
this figure the vertical edges are enhanced and the horizontal edges are also 
preserved. 


2. Recursive Implementations; The basic difference between 
nonrecursive and recursive implementations is that in recursive implementations 
the output computed at any stage is dependent on at least one of the previously 
computed values of the output. That is, the equation giving the output is of 
the form 


M-l 

£ h j 

j~o 


L-l 

y i-j - £ § k *i. k 

k=0 


(3.38) 


where hg = 1 (with no loss of generality) and there exists at least one j 
such that 1 < j < M-l and hj - 0. 

It has been shown by Rabiner and Schafer [l8 ] that narrow band filters 
designed by the frequency sampling techniques (section III.A.l.e) can be 
implemented recursively faster than the fast convolution method. 


In this section we describe another class of filters which lend them- 
selves to very fast recursive implementations even in two dimensions and 
also show how they can be used to generate matched filters. 

a. Spatial and frequency domain characterization of the 
filters: First of all, consider the one -dimensional filter which replaces 

every pixel density by the average of the densities of 2L + 1 neighboring 


pixels including the pixel itself.' That 
and yi is the output sequence, then 



/ L 


Yi - 

(k?L X i+k 

) /(2L+ 

Thus, the filter 

weights are 

gi where 

gi 

1/(2L + 1) 

for 

= 

0 

for 


is, if x^ is the input sequence 
) (3.39) 

i€ [ -L, L] 

if [-L, L] (3.40) 
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Figure 3,16 Band Emphasis Filtered Test Pattern 
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This filter is a linear phase filter whose modulation transfer function 
is shown in Figure 3.17 for L=3. It can be seen from Figure 3.17 that the 
'major lobe' is at the low frequency end and the first zero is at l/(2Lfl) 
times the sampling frequency. Thus the filter is, in effect, a low pass 
filter, even though it is not a good approximation to an ideal low pass 
filter. 


Now, from the above low pass filter it is easy to obtain a high 
pass filter by subtracting the output of the low pass filter from the 
input. That is, the average density of (2L + 1) neighboring pixels is 


subtracted from each pixel density, 
filter, we can write 

If 

hi are 

the weights of this 

h t = <5i - 1/ (2L 4- 1) 

for 

ie [-L, 

L] 

= 0 

for 

i i [-L, 

l] 

where 




(5i = 1 for i = 0 and 

*i" 

0 for i 

4 o. 


(3.41) 


The modulation transfer function of this filter is shown in Figure 3.18 for 

L = 20. 


In general, the filter transfer function is derived from the filter 
weights by employing the discrete Fourier transforms such that the (complex) 
transfer function is 

L 

G (^k) = E g ± exp(-jidk u/ 0 ) (3.42) 

i= "*Ii 

where 4 is the sampling interval in the original domain and u>q is the sampl- 
ing interval in the transform (frequency) domain. The zero frequency 
behavior of the filter is given by G(wk) evaluated for k = 0. Ideally, a 
low pass filter has unit transmission at zero frequency while the high 
pass filter has zero transmission at this frequency. Thus general constraints 
on the filter weights are 

L 

for a low pass filter (3.43) 


for a high pass filter (3.44) 

and these constraints are seen to be satisfied by the filter weights of 
Equations (3.40) and (3.41). 


a. 

i=-L 

L 

E 

i=-L 


= 1 


Si = 0 
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FIG1. LP FILTER WITH L=3 



Figure 3.17 Low Pass Filter with L=3 
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FIG2. HP FILTER WITH L = 20 
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Figure 3.18 


High Pass Filter with L=20 
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In discussing the design of the class of filters based upon constant 
weight low pass and high pass elements, it is convenient to represent the 
filter function graphically in terms of the signs of the weights. This 
representation is shown in Figure 3.19 for the basic low pass and high pass 
filter functions, and it will be employed extensively in the remainder of 
this discussion. Filter sample values may be assumed to be zero if no sign 
is associated with the sample. It should be obvious from Figure 3.19 that 
similarity of sign does not imply equality of weight, but the value of weight 
to be associated with a given sign can generally be evaluated simply. 


The basic low pass filter of Equation (3.40) can be used to produce 
a band pass filter using two approaches. One approach is to have two low pass 
filters g! and gV as in Equation (3.40) of lengths (2K + 1) and (2L + 1) 
respectively where K<L and to subtract the output of the second from that 
of the first. This, in effect, produces a filter whose weights are 


h i ." 


2K + 1 


2L + 1 


for i € [-K, K] 


for i e [-L, L ] - [-K, k] 
for i 4 [-L, L ] 

The modulation transfer function corresponding to the weights in 
Equation (3.45) is shown in Figure 3.20 for K = 3 and L = 20. 


= _ 1 

2L + 1 

= 0 


(3.45) 


Another approach for generating a band pass filter is to use a low 
pass filter (or Equation (3.40)) and a high pass filter (of Equation (3.41)) 
in tandem. Then, the modulation transfer function of the band pass filter 
will be the product of those of the low and high pass filters. The filter 
weights in this case are obtained by 


00 

“ ]C §k h 'i-k for a11 1 (3.46) 

k= - oo 

where | g^J are the low pass filter weights given by (3.40) with K instead of 
L and|hiP are the high pass filter weights of Equation (3.41) and K< L. It 
is easy to see that 


h i 


1 

2K + 1 
1 

2L + 1 
1 

2L + 1 
0 


1 — for i e [ -K, K 1 

2L + 1 L ’ J 

for ief-K-1, -L+k]u[k+1, L-k] 

1 - . jl*] . -< L ~ K ) ) for if[-L-K, -L+K- 1 ] y [L-K+ 1 , L+K ] 

for i 4 [ -L-K, L+K] (3.47) 


The modulation transfer function of this type of band pass filter is 
shown in Figure 3.21 for K = 3 and L = 20. 
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FILTER FUNCTION FILTER FUNCTION 

+ + + + + + + + + + + + + _ ^ _ + _ __ __ _ 
REPRESENTATION BY SIGN OF WEIGHTS REPRESENTATION BY SIGN OF WEIGHTS 

(a) LOWPASS FILTER (b) HIGH PASS FILTER 

FIGURE 3.19 REPRESENTATION OF CONSTANT WEIGHT FILTER FUNCTIONS 
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Figure 3,20 Band Pass Filter as the Difference between two 
Low Pass Filters 
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FIG4. BP FILTER WITH LP AND HP IN TANDEM 
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Figure 3.21 Band Pass Filter with Low Pass and High Pass 
Filters in Tandem 
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Filter realizations: 


b . 


(1) Z - Transform representation of filter functions: A 

convenient way to obtain the recursive realization of a filter in general 
is to examine the Z-transform of its impulse response. If the impulse 
response (represented numerically by the filter weights) of a filter is 
given by the sequence jgi}, then the Z-transform G(Z) is given by (see [19], 
for example) 

00 

G(Z) = 23 g-jZ" 1 (3.48) 

i=— 00 

which is defined for all Z if jgi} has only a finite number of nonzero ele- 
ments . 


Now, the Z-transform of the output {y^j is related to the Z-transform 
of the input {x^} by 


Y(Z) = G(Z) X (Z) 


(3.49) 


Therefore, from Equation (3.40), the input/output relationship for the low 
pass filter can be written as 


Y(Z) 


i (z L - z- tt* 1 )) 

(2L + 1) (1 - Z-l) 


X(Z) 


(3.50) 


The realization of this as a first order recursive filter is shown in Figure 
3.22. Here Z“1 denotes unit delay. The Z-transform of the impulse response 
of the high pass filter of Equation (3.41) is given by 

1 (Z L - Z’ (L+1 >) 

(2L + 1) (1 - Z" 1 ) (3.51) 


H(Z) = 1 - 


The recursive realization corresponding to Equation (3.51) is shown in Figure 
3.23. Here, LP(L) denotes the low pass filter shown in Figure 3.22. 

It is easy to see that the band pass filters defined by Equations 
(3.45) and (3.47) can be realized as shown in Figures 3.24 and 3.25 respectively. 
Also, recursive realizations can be readily derived for some two dimensional 
low, high and band pass filters from the basic low pass filter realization 
in Figure 3.22. 


(2) Two dimensional filters: Two general types of two-, 

dimensional filters are realizable by a 'two pass’ procedure, as discussed 
in the beginning of Section III-B: (i) "Sum" filters whose transfer func- 

tion can be written as a sum of two functions each being a function of a 
single frequency variable; (ii) "Product" filter whose transfer function can 
be written as a product of two functions each being a function of a single 
frequency variable. "Sum" filters can be implemented by filtering the input 
separately in the horizontal and vertical directions and then adding the 
outputs. "Product" filters can be implemented by filtering the input in one 
direction and then filtering the result in the perpendicular direction. 

Some of these two-dimensional filters are discussed below. 
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FIGURE 3.22 RECURSIVE REALIZATION OF THE LOWPASS FILTER OF EQUATION (3.40) 



FIGURE 3.23 RECURSIVE REALIZATION OF THE HIGHPASS FILTER OF EQUATION (3.41) 
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FIGURE 3.24 RECURSIVE REALIZATION OF THE BANDPASS FILTER OF EQUATION (3.45) 


61 










FIG. 3.25 RECURSIVE REALIZATION OF THE BAND PASS FILTER OF EQUATION (3.47) 



FIG. 3.26 RECURSIVE REALIZATION OF A 2 DIMENSIONAL LOW PASS FILTER (MOVING 
AVERAGE GENERATOR ) 



FIG. 3.27 ALTERNATE REPRESENTATION OF 2 DIMENSIONAL LOW PASS FILTER 
IN FIGURE 3.26 
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Suppose, for example, that it is necessary to implement a filter 
whose output corresponding to a given point in a picture is the average 
density of a (2L-|+1) x (2L2$-1) rectangle centered about the given point. 
Then, denoting basic low pass filters in the horizontal and vertical 
directions by LPjO and LJ^C.) respectively, the above filter can be 
realized as shown in Figure 3.26. Of course, there is a difference in 
the operations of LP^(.) and LP 2 (.) in that LP-^(.) handles points within 
a record (it operates along rows of the picture array) while LP 2 (.) 
handles points from several records (it operates down columns of the 
picture array). Table 3.4 shows the realizations of a number of two- 
dimensional filters using the basic low pass filter of Figure 3.22. The 
table also shows the nature of the point spread function of each filter, 
represented by the signs of the weights. In the implementation of the 
filter shown in Figure 3.26, for example, the tape-to-core and core-to- 
tape transfers and tape rewinds can be minimized by having the data to be 
filtered on two tapes and reading the appropriate records from each tape 
and then operating LP^O.) on them. These operations are represented in 
Figure 3,27. Here Z 2 “l is a unit delay in the vertical direction. Thus, 
when the i fc h record is being filtered, the (i + l^)*-^ and (i - L 2 -l) st 
records are read, each from a different tape, and their difference is 
added to the (i - l) st record of y' which is held in core. LPj(Li) is 
then operated on y'. 

(3) Practical aspects of filter implementation: The manner 

of implementing the recursive realization of the filters summarized in Table 
3.4 is dependent entirely on the medium available for storage of the image 
data. If sufficient core storage is available to store the entire image 
array, the entire filter functions may be implemented as algorithms. 

More realistically, the image array is stored on drum, disk or tape, 
and to avoid excessive data seeking, the algorithms are modified slightly 
and the image data records are restructured to facilitate processing. 

The Z-transform variable may be identified with the forward shift 
operator E associated with finite-difference equations. A filter operator 
Z^ may thus be interpreted as a forward shift of a data sequence by L sample 
intervals, and Z“ L may similarly be interpreted as a backward shift by L 
sample intervals (or correspondingly a delay of L units). When a data 
sequence is stored in its entirety in core, a forward shift with respect 
to a specific sample may be realized almost instantaneously, but this is 
not feasible when the required data sequence is held within records stored 
sequentially on tape. The forward shift operation is then implemented 
physically by duplicating the data records on tape and advancing the 
supplementary tape by the appropriate shift interval. 

It is found in general that the number of arithmetic operations is 
reduced if the vertical filtering is done first. The number of tape rewinds 
is minimized by having the data on several tape units. The number of tape 
units on which input data is to be present depends on the number of input 
records needed simultaneously to form one record of output. 
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TABLE 3,4 RECURSIVE FILTER REALIZATIONS 









































x I a 



+ 

+ +I+-*- 

+■ 

+ + 1 + ■t- 


+ + 1 + -*- 

1 1 + II 

1 1 + 1 1 

+ 

+ +■ 1 + + 

+ 

+ + 1 + + 

+ 

+ + 1 + + 




65 


TABLE 3.4 RECURSIVE FILTER REALIZATIONS 
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TABLE 3.4 RECURSIVE FILTER REALIZATIONS 


















































Also, the edges of the picture should be augmented properly in order 
to avoid unpleasant edge effects. This can be done by defining the data 
outside the domain of the picture arbitrarily (See Section III-B.l). It 
has been found from previous experience that good visual effects are obtained 
by defining the data outside the domain to be the reflection of the data in 
the nearest edge. The reflections of data needed while filtering each record 
(i.e., horizontal filtering) can be generated in core. But the reflections 
needed for vertical filtering have to be stored on tape. The number of 
records which should be reflected at the top (and bottom) of the picture is 
K where (2K + 1) is the maximum anticipated vertical filter length. In fact, 
it is convenient to handle the records if a tape is prepared with the 
reflection at the top appearing first, the data records themselves appearing 
next and the reflection at the bottom appearing at the end. 

In Table 3.4, the number of tape units on which input data should be 
supplied and also the number of tape units on which data and reflections 
should be supplied are shown for each of the filters. 

c. Extension to matched filter implementation: We shall now 

present a simple extension of the recursive filters discussed above to some 
particular cases of matched filters. We shall restrict ourselves to cases 
where only the variations of density over the local average are of interest. 

We can therefore choose the weights g^£in equation (3.26) such that 

E 0 0-52) 

k, l 


(this condition is satisfied by filters which are not low pass at least in 
one dimension). Also, equation (3.26) is now modified to 


where 



Xi J (2K+ 1) (2L+ 1) ^ i X i+M+£ 

Now, from (3.52) it follows that the numerator of (3.53) is 


(3.53) 


(3.54) 


K L 

H J2 s k£ x i+k, j+£ 

k=-K £=-L 

Also, from (3.54), it is easy to see that 

2 

E (x i+k,j+£- ; ij> ■ E ’‘l+k.Jfi- < 2K «> (3-55) 

k, & k , L 
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Further, note that §k £ is a constant over all i,j and hence 

if linear scaling will be used for displaying data (as is most often the 
case) this factor can be omitted from the expression (3.53) for j • 

Therefore, the recursive filters described above can be used to 
generate y.. with the modification shown in Figure 3,28. Note that in all 
the filters' 1 shown in Table 3.4, U^O^) appears either explicitly or 
implicitly and that in some cases the tandem combination I^CL^LP^CL^) on 
x in Figure 3.28 can be eliminated. Also, no additional input/output 
operations are involved in the computation of the matched filter output 
given by (3.53) compared with the computation of only the numerator of (3.53) 
as is done in the case of the filters in table 3.4. Thus we have an 
efficient implementation of matched filters which can be used for detection 
of shapes resembling the point spread functions in table 3.4. 

Examples of matched filtering of a test pattern using the above method 
are shown in figures 3.29, 3.30 and 3.31. The basic test pattern used in 
these figures consists of vertical lines of high density (density number = 15) 
separated by areas of low density (density number = 0). The vertical lines 
have widths 1, 5, 9, ..., 25 pixels and the areas separating them have widths 
of 25 pixels. The template used in all these pictures is a 39 x 49 pixel 
rectangle with a vertical line of width 13 pixels at the center. 

In figure 3.29, the frame in the top left comer shows the test pattern 
rescaled between 0 x 63. The second frame in the first column is the result 
of matched filtering which is linearly rescaled so that the minimum is set 
to 0 and the maximum is set to 63. The remaining frames are all results of 
rescaling the filtered output with progressively increasing thresholds, by 
setting values below the threshold to 0 and linearly scaling the values 
above the threshold to lie between 0 and 63. 

In figure 3.30, the frame in the top left comer shows the unfiltered 
test pattern corrupted with additive noise which is approximately white and 
Gaussian with a mean of 15 and variance of 8. The noisy pattern is rescaled 
linearly to lie between 0 and 63. The second frame and all the other frames 
are the results of matched filtering, thresholding and rescaling as in 
figure 3.29. 

Note that as the threshold is increased in figures 3.28 and 3.29, the 
vertical lines get thinner and thinner, the lines farther away from the 
fourth vertical line (i.e. the line of width 13 pixels) disappear first and 
in the last (bottom righthand corner) frame only the line corresponding to 
the fourth line remains. This clearly demonstrates the effect of mismatch 
upon the matched filter output. 

Figure 3.30 shows the effect of varying the input noise level on the 
matched filtered output. The left column shows the test pattern corrupted 
with increasing noise levels as one proceeds downward. The right column 
shows the corresponding matched filtered outputs which have been thresholded 
and rescaled. The threshold is taken as 0. That is, all output values 
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FIGURE 3.28 MATCHED FILTER REALIZATION 
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Figure 3,29 Test Pattern and Matched Filter Output with 
Different Thresholds 
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Figure 3.30 Noisy Test Pattern and Matched Filter Output 
With Different Thresholds 
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less than 0 have been set to 0 and the nonnegative values have been linearly 
rescaled to lie between 0 and 63. It can be noted that the filtered pictures 
are reasonably good replicas of the original (noise free) picture except that 
in the high noise cases (last two) there is some noticeable distortion. 

f. Comparison of implementation times: It might be of some 

interest, to compare the times taken to filter a picture using nonrecursive 
implementation using fast convolution and the recursive implementation of 
the particular filters discussed above. Of course, it is to be remembered 
that the recursive implementation discussed in Section III-B.2 is only for 
particular transfer functions and the nonrecursive implementation in Section 
III-B.l can implement any transfer function. Table 3.5 shows the typical 
values of the time needed to filter a 256 x 256 picture in one dimension 
using different lengths of filters. 


Filter Length 

Time (Min.) 

40 

15 

150 

40 

200 

84 

240 

85 


Table 3.5 Times for Implementation of a 
General One -Dimensional Filter 
(Nonrecurs ive) 

These times can be reduced slightly if the fact that the data and the filter 
weights are real is exploited. However, their orders of magnitude will be 
retained. Further, if two-dimensional filtering is needed, a two pass pro- 
cedure has to be employed (under the present core limitations) involving 
two large matrix transpositions and two one -dimensional filtering operations. 
Consequently, the times shown in Table 3.5 will be more than doubled for 
two dimensional filtering. 

In the case of radiographic images in which fine features such as 
minute cracks are to be detected it is necessary to digitize the imagery at 
the highest avilable resolution, viz., 12. 5 fl sampling interval. This means 
that even for an image size of 12.5mm x 12.5mm a 1000 x 1000 digital array 
is generated. Therefore, the nonrecursive filters mentioned above become 
rather impractical. 

Experiments performed on a radiographic image of size 1000 x 320 and 
a filter which was band pass in the horizontal direction and low pass in the 
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vertical direction (see table 3.4) with parameters kx = 6, Lx = 24 and L£ = 19 
have shown that the time taken for recursive implementation in two dimensions 
was of the order of 8.5 minutes which is a considerable reduction from the 
times shown in table 3.5. With the same filter and picture sizes the matched 
filter implementation took about 14.2 minutes. 

3. Sequential similarity detection |2o] : A promising new technique 

for the "automatic determination of the local similarity between two structured 
data sets" has been reported by Barnea and Silverman. They also illustrate 
the application of this technique for registration of remotely sensed imagery. 
This technique, called Sequential Similarity Detection works on the principle 
that similarity between an m x n template and an area on a picture can be 
defined in terms of an error criterion which is a sum of mn positive terms, 
each of the terms being the error between the template and the picture area 
at a particular point. Thus, the error is a monotonic nondecreasing function 
of the number of points at which it has been evaluated and summed. It 
increases slowly when the template matches the picture area approximately 
and increases rapidly if the template does not match the picture area. Thus, 
if a threshold is properly chosen on the error, one can decide whether a 
match occurs or not without actually computing the error completely for all 
the points. This is particularly advantageous in cases where the number of 
points at which the match occurs is small compared to the total number of 
points in the picture. Barnea and Silverman claim that this method is about 
50 times faster than even FFT methods of correlation for registration. 

The applications of this technique to image enhancement are being 
tried at present. Experiments are being performed on test patterns to 
develop criteria for threshold selection so that the processing time is 
minimized and the probability of detection of desired features in a corrupted 
image is maximized. The results of these experiments will be reported in 
subsequent reports. 
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SECTION IV. APPLICATIONS 


The image processing techniques and software discussed in the 
previous sections have been applied to several radiographs. It is the 
purpose of this section to demonstrate these applications. Discussed are: 
(i) the enhancement of flaws in the radiograph of a structural weld and (ii) 
the enhancement of the radiograph of a metal with a crack in it in order to 
determine the extent and shape of the crack. 


A. Structural Weld Radiograph 

Figure 4.1 shows the radiographic image of a structural weld. This 
image was prepared by first digitizing the radiograph with 25 H scanning 
interval and then writing it on film using the Photowrite and printing the 
resulting image. While there are three flaws which are clearly visible 
(one large approximately circular dark spot, one approximately elliptic 
dark spot and a small spot about 1/2" to the right of it) it is not clear 
whether there are any cracks joining these flaws. The purpose of enhance- 
ment was to make the shapes of the flaws more evident and to detect cracks, 
if any, between the flaws. 

A band emphasis filter was applied to the image in figure 4.1 so 
that edges and cracks would be enhanced but the high frequency noise would 
be suppressed. The transfer function used was the one shown in figure 3.6 
except for the fact that the gain in the attenuated bands was set to .1 
instead of 0. Figure 4.2 shows the effect of one dimensional filtering of 
the region containing the flaws in figure 4.1. The filtering was horizontal 
so that the vertical (and almost vertical) edges of the flaws have been 
enhanced. Figure 4.3 shows the effect of two dimensional filtering of the 
same region. The filter was implemented as a two pass procedure, as illustrat- 
ed in figure 3.7. The transfer function was the same band emphasis function 
as above in both the horizontal and vertical directions. In this case both 
the horizontal and vertical edges are enhanced. Note the difference between 
the edges of the two larger flaws in figures 4.2 and 4.3 and in particular 
the fact that the nearly horizontal edge of the oval flaw is marked much 
better in figure 4.3. 

The above filtered images show the shapes of the flaws quite clearly 
and indicate no crack joining the two larger flaws. However, there is some 
doubt as to whether the oval flaw and the small flaw to its right are joined 
by a hairline crack. In order to confirm the existence (or otherwise) of a 
crack, the region including the oval flaw and the small flaw was digitized 
with 12.5 scanning interval. The scanned area is shown in figure 4.4. 

The scanning direction in this case is perpendicular to that in figure 4.1 
and the oval flaw appears vertically above the small flaw (instead of on 
the left). Figure 4.5 shows the effect of contrast enhancement by point 
operations of linear density stretching. (see section III-A.2.b). A 
256 x 256 pixel area of figure 4.4 containing the flaws of interest was 
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Figure 4.1 Radiographic Image of a Structural Weld 
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Figure 4.2 Image after Horizontal Band Emphasis 
Filtering of Part of Figure 4.1 
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Figure 4.3 Image after Two Dimensional Band Emphasis 
Filtering of Part of Figure 4.1 
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Figure 4.4 Small Part of Figure 4.1 Digitized with 
12.5/1 Scanning Interval 
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Figure 4.5 


Image after Linear Density Stretching was Applied 
to Part of Figure 4.4 for Contrast Enhancement 
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extracted and the minimum density number in the area was subtracted from 
the density number at each point and the resulting numbers were linearly 
rescaled between 0 and 63. Enlargement by a factor of 3 in each direction 
was obtained by repeating each pixel density three times in both the hori- 
zontal and vertical directions. The flaws are very clearly visible in 
figure 4.5 and this figure indicates that there is no hairline crack join- 
ing the two flaws. This area was further explored by one dimensional and 
two dimensional filters having several different transfer functions. 

Figure 4.6 shows six cases of one dimensionally high pass filtered images. 
The filtering in all the cases was in the horizontal direction. The 
implementation of the filters was recursive as discussed in section III-B.2 
and the lengths of the filters for the six cases are 3, 23, 43, 63, 83 and 
103 pixels respectively. The filtered data was linearly rescaled between 
0 and 63 in all cases. The first three cases are shown in the first column 
and the last three in the second column of figure 4.6. The filter length 
increases with each frame as one proceeds downward in the figure. It can 
be seen that the smaller the filter length, the sharper are the vertical 
edges of the flaws except in the case of filter length 3 (top left frame) 
where the high frequency noise is predominant. Figure 4.7 shows the effect 
of two diminsional high pass filters of various sizes on the same region as 
in figures 4.5 and 4.6. A 1 x 1 'filter' was also attempted just as check 
on the program. This corresponds to subtracting the density at each pixel 
from itself and hence, as expected, the top left corner of figure 4.7 con- 
sists of a 'frame' with all densities equal to 0. The other five frames 
correspond to filter sizes of 11 x 11, 21 x 21, 31 x 31, 41 x 41 and 51 x 51 
respectively, the order of frames being the same as in figure 4.6. As in 
figure 4.6, the shorter filters yield sharper edges and in this case the 
comment applies to both horizontal and vertical edges of the flaws. From 
these pictures it can be seen that there are no hairline cracks between the 
flaws . 

B. Detection of a Crack in a Metal 
from its Radiograph 

Figure 4.8 shows the radiographic image of a metal. This image was 
obtained by scanning the radiograph with a 12.5// interval and writing it 
back using the Photowrite and printing the image. The lighter rectangular 
area is the region of interest. A thin vertical crack can be seen to begin 
at the top of the picture. The crack is visible better in the radiograph 
than in the printed version of it. However, the end of the crack cannot 
be seen either in the original radiograph or the printed version. The pur- 
pose of enhancement was to find the extent of the crack by making the crack 
more prominent with respect to the background. 

Figure 4.9 shows the result of subtracting the minimum density in 
figure 4.8 and stretching the densities between 0 and 63. It can be seen 
from this figure that the crack (the dark vertical line at the top of the 
picture) has become more prominent and merges into a dark patch on the right 
half of the picture (about 2.5" from the top) and extends beyond the patch 
sloping towards the right. 
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Figure 4.7 Effect of Various Sizes of High Pass 

Filter (Two Dimensional) on Figure 4.4 
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Figure 4.9 


Linearly Rescaled (Density 
Version of Figure 4.8 


Stretched) 
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Figure 4.10 shows the result of applying a matched filter which en- 
hances vertical lines to figure 4.9. The filter (described in section 
III-A.2.a) gives a high output for the positions in the picture which 
resemble a template. The template used in this case was a 39 x 49 pixel 
rectangle consisting of a 13 pixel wide vertical line in the center, so 
that the filter enhances vertical lines darker with respect to the back- 
ground. A comparision of figures 4.9 and 4.10 shows that the crack in 
figure 4.9 yields a very dark vertical line (the darkest line beginning 
at the top of the picture). 

Figure 4.11 is the result of thresholding the matched filter output 
before rescaling it between 0 and 63 so that only the regions of the picture 
which are, on the average, darker than the background along a vertical line 
are retained and the other points are set to 0 density. The position of the 
crack and the direction in which it proceeds are clearer in figure 4.11 than 
it is in figures 4.9 and 4.10. Figures 4.12, 4.13, 4.14 and 4.15 demonstrate 
the same set of operations performed on another radiograph. The crack is 
again more prominent and clearly visible in the filtered images. Figure 4.16 
shows the effect of adding the densities in the rescaled version of figure 
4.13 and twice the densities in the filtered and thresholded version of figure 
4.15, point by point. This picture consists only of approximately 5/9 t ^ ls from 
the top of the pictures in figures 4.12 through 4.15 and has been enlarged. 

The crack is more prominent in figure 4.16 than in figure 4.13 and the back- 
ground which was missing in figure 4.15 has been restored. 
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Result of Filtering Figure 4.9 with a Matched 
Filter to Enhance Dark Vertical Lines 
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Figure 4.11 Result of Thresholding the Matched Filtered 
Output in Figure 4.10 Before Rescaling 
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Figure 4.13 Linearly Rescaled Version of Figure 4.12 
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Figure 4. 14 


Result of Filtering Figure 4.13 with a Matched 
Filter to Enhance Dark Vertical Lines 
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Figure 4„ 15 Result of Thresholding the Matched Filter 
Output in Figure 4„ 14 Before Rescaling 
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Figure 4.16 Result of Point-by-Point Addition of 
Densities in Figures 4.13 and 4.15 
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SECTION V CONCLUSIONS 


The application of existing techniques of digital signal processing 
to radiographic image enhancement has been considered. A versatile set of 
programs has been developed to perform the image handling and image process- 
ing operations. These programs have been designed on the basis of a tape- 
drive-based IBM 7094. 

The main image handling operations considered are: (i) converting the 

digitized data generated by a microdensitometer into a format suitable for 
use on the computer, (ii) automatic extraction of the data corresponding to 
the region of interest from a tape generated by the microdensitometer, (iii) 
transposition and 90° rotations of large data arrays (which do not fit in 
core) (iv) translation of data arrays for registration (v) reduction of the 
dimensions of data arrays by integral factors and (vi) reconverting the data 
from computer format to a format suitable for "Photowriting" or display so 
that several frames of data can be centered and written on film or displayed 
on a screen. 

Both the frequency domain and the spatial domain approaches have been 
considered in designing the image processing software. Programs have been 
developed using the Fast Fourier Transform for the design of filters with 
arbitrary continuous transfer functions to within a specified error, so 
that the .lumber of filter weights is as small as possible. These programs 
are useful for designing digital filters which can enhance, attenuate or 
stop any desired ranges of spatial frequencies in the pictures. Fast imple- 
mentation of the digital filters have been developed using fast convolution 
techniques and for some special cases of digital filters using recursive 
techniques. It has been found that the recursive implementation is much 
faster in the cases of interest in the present work, particularly for 
enhancing vertical or horizontal lines in radiographs. In these cases 
matched filters (template matching filters) have been implemented very fast 
using the recursive method (taking about 14 minutes for filtering a 1000 x 
320 array) . 

Examples of test patterns have been used in this report to illustrate 
the techniques and the application of the techniques to radiographic images 
has been shown through several examples. 

The choice of filter characteristics -transfer functions in the 
frequency domain, template sizes and certain thresholds on normalized cross 
correlation in the case of matched filtering - is still empirical. The 
way they are chosen at present is to filter a small area of the image with 
several filters and choose the filter which gives the best (visual) results. 

A knowledge of the statistical characteristics of noise in the radiographic 
imaging systems will help in the choice of the filter parameters. But the 
results obtained by the empirical approach are quite encouraging. 
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It is felt that further developmental effort should be in the areas 
of template matching and other particular-feature-oriented approaches. Cases 
in point are a fast approach for registration using sequential similarity 
detection [ 20] which can be applied to template matching and image pattern 
recognition and "Edge and Curve Detection" reported by Rosenfeld et al [ 2l] 
(July 1972) which are used for the detection of "texture edges" in images. 
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PART II 


SECTION I. I-NTRODUCTION 


This part of the report presents the formal documentation of the 
image handling and processing software discussed in Part I. Section II 
presents the details of the programs for image handling operations such 
as automatic picture extraction, transposition and 90° rotations of pic- 
ture data, generation of reflections of picture data about the top and 
bottom edges arid packing of data for multiple image display. Section III 
presents the details of the programs for image processing operations 
including design and implementation of filters in frequency domain, point 
operations on pictures such as rescaling and density manipulation and 
fast recursive implementation of a particular class of filters. Complete 
optimization of the programs has not yet been attempted and some of the 
programs can be improved considerably. 


SECTION II. IMAGE HANDLING PROGRAMS 


Al. Automatic Picture Extraction 


1. NAME 

Deck Name: APEXS 

2. PURPOSE 

To determine the region of interest in a digitized picture, extract 
it, unpack it and linearly rescale the data (as an option) and write it 
in unpacked format on a tape. Also, as an option, to override the 
determination of the region of interest and just extract, unpack and 
write a specified rectangular area of a digitized picture. 

3. CALLING SEQUENCE 

CALL APES (IA, IX, ILIM, IAMX, IAMN) 

where all the calling arguments are integer work arrays to be dimensioned 
as described in the comments in the listing attached. 

See also paragraph 7.3 below. 
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INPUT -OUTPUT 


4.1 Input: The input tape should be on logical unit 8. (See note 

on buffer size for unit 8 in the comments in the listing.) The input 
tape is the tape produced by the microdensitometer and is assumed to 
contain integers between 0 and 63 packed as 6 pixels per computer word. 

4.2 Output: The output of the program will be on logical unit 10, 

written in FORTRAN binary format, each record consisting of the record 
number followed by the data in unpacked format. 

4.3 File Storage: None. 

5 . EXITS 


When the .program is used in the mode where the region of the picture 

to be extracted is to be determined, the program exits giving the mes- 

sage "NO RECORDS WERE EXTRACTED SINCE THE DENSITIES WERE ALL .GT. 62" 
if either there are no records in the input file with density numbers 
less than 63 or the number of successive occurrences of such records 
after the first such record has been found is less than or equal to the 
vertical margin MVl (see paragraph 7.3 below). 

While determining the region to be extracted, the program finds the 
first and last points to be extracted in each record. The initial column 

number is taken to be the right -most of the initial points of all the 

records to be extracted and final column number is taken to be the left- 
most of the final points of all the records to be extracted (see tech- 
nique description in Part I, Section II -A, for details). If the initial 
column number exceeds the final column number (after margins are allowed), 
then the program exits, printing the message "NO RECORDS WERE EXTRACTED 
SINCE INITIAL COLUMN EXCEEDED FINAL COLUMN." In this case, the "frame 
skeleton details," i.e. , the initial and final points in each record 
which are regarded by the program as belonging to the region of interest 
(after allowing specified margins) will be printed also. 

In either of the above cases, no records will be written on unit 10. 
The messages will be printed both on-line and off-line. The frame 
skeleton details will be printed only off-line. 

6. USAGE 

The program is written in FORTRAN IV. It is presently implemented 
on IBM 7094. 

7. EXTERNAL INTERFACES 

7.1 System Subroutines: SYSLOC, SKRBIN, BSFILE 

7.2 Other Programs Called: UPVECT, UPVEC and BECOL (See Sections 

II. B and II-C.) 
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7.3 


External Storage Used: 

(1) COMMON/ IMPE/HSL, SINT ,MH1 ,MH2 ,MVl ,MV2 , ISCAL, IFLG 

All parameters in this common block are inputs. The 
definitions of each parameter are given in comments in the listing. 

(2) COMMON/ OVRID/lRINIT , IRFIN, ICI, ICF , IFLG1 

All parameters in this common block are inputs. The defi- 
nitions of each parameter are given in comments in the listing. 

Note: If the determination of the region of interest is to 

be overridden (i.e., IFLGl^O) , only HSL and SINT need 
be specified in the common block /INAPE/. If IFI£1=0, 
then no parameter other than IFLGl need be specified 
in the common block /OVRID/. 

(3) COMMON/ OUAPE/ NREC , NEL 

Both the parameters are outputs giving the number of 
records and number of pixels per record which were actually extracted 
and written on unit 10. 

8. PERFORMANCE SPECIFICATIONS 

8.1 Storage: 2176g locations of core. 

8.2 Execution time: The execution time is highly dependent on the 

size of the image scanned, and the size of the window being extracted. 

It also depends on the location of the window refative to the first 
record. The following figures are quoted as an example. In the test 
case attached, it can be seen that the time for determining the region of 
interest, unpacking and rescaling was about 7.5 minutes including tape 
rewinding at the end of the job. 

8.3 I/O load: Under conditions not covered in paragraph 5 above, 

the program produces an output tape on logical unit 10. The values of 
the initial row, initial column, final row and final column extracted 
will be printed both on-line and off-line. The histogram of the den- 
sities in the region extracted will be printed both on-line and off-line. 
If scaling is requested (i.e., ISCALfO) , the histograms of both the 
original and rescaled versions will be given in the off-line output. 

If IFLG^O, then the "frame skeleton" details and the minimum and maxi- 
mum value of the densities in each record in the region of interest will 
be printed. 


8.4 Restrictions: None. 
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9. 


METHOD 


See Part I, Section II-A, for a description of the method. 

10. COMMENTS 

The program might give the second message discussed under paragraph 5 
due to skewness of the rectangular region of interest with respect to the 
scanning direction. In this case, it is suggested that the program be 
retried with either manual override (IFLGl^O) or with increased vertical 
margins (MVl and MV2). 

11. LISTING 

A listing of the program is attached at the end of Section II-A3. 

12. TESTS 

— i 

The program has been tested on many digitized data tapes which were 
generated by masking all but the desired rectangular regions in several 
radiographs. It hks been found to work satisfactorily in all cases where 
most of the density numbers were less than 63. Typical values of margins 
to be used are 25 to 40 pixels, the required margins being larger for 
smaller scanning intervals. A typical output of a successful run is 
shown at the end of the listing. 

I 

\ 
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A2. Unpacking a Specified Part of a Record 


1. NAME 

Deck name : UNPVEC 

Entries: UPVEGT and UPVEC 

2. PURPOSE 

To read a record of packed data and unpack a specified part of it 
(UPVECT) or to unpack a specified part of a given record (UPVEC). 

3. CALLING SEQUENCE 

CALL UPVECT (INITC, NEL, IX, IA, N, NERR) 
or 


CALL UPVEC (INITC, NEL, IX, IA, N) 

IX is the integer work array into which the packed record is read. 

IA is the integer output array into which the desired part of IX 
will be unpacked. 

INITC = input integer giving the unpacked word number in IX which 
should go into IA(1). (See paragraph 9 below.) 

NEL = input integer giving the number of (packed) words to be read 
into IX from the input tape (UPVECT) or the number of words in IX (UPVEC). 

N = input integer giving the number of words desired in the output 
array IA. 

NERR = output integer equal to 0 if the routine does not read an 
end of file and equal to 1 if it does. 

4. INPUT-OUTPUT 

4.1 Input: The input data for entry UPVECT is assumed to be on 

logical unit 8. The buffer size for unit 8 should be at least NEL. The 
tape is assumed to contain integers between 0 and 63 packed as 6 pixels 
per computer word. For entry UPVEC, the input is through the calling 
sequence only. 

4.2 Output: The output is through the calling sequence only. 

4.3 File storage: None. 
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EXITS 


When the program UPVECT encounters an end of file mark, it sets all 
words in IA equal to 0, prints the message "EOF ENCOUNTERED BY UPVECT" 
both on-line and off-line and returns. 

6 . USAGE 


The program is written in FORTRAN TV. It is presently implemented 
on IBM 7094. 

7. EXTERNAL INTERFACES 

7.1 System subroutines: SYSLOC, REDTPR 

7.2 Other programs called: UNPINT (a program which unpacks a given 

word into an array of 6 integer words) . 

7.3 External storage used: None. 

8. PERFORMANCE SPECIFICATIONS 

8.1 Storage: 354g core locations 

8.2 Execution time: TBD 

8.3 I/O load: None. 

8.4 Restrictions: None. 

9. METHOD 

The program UPVECT reads NEL words of a record of data from the 
input tape into IX, using REDTPR. When unpacked, the array IX would thus 
yield 6 NEL words. It is desired to transfer N words out of these 6 NEL 
words into the array IA, starting from the INITCth word. The program 
determines the word number in IX corresponding to the INITCth unpacked 
word and unpacks only the number of words required to produce N unpacked 
words. Also, if (N + INITC - 1) exceeds the number of words NW read 
from tape (as determined by the routine REDTPR) the last few words of IA 
will be set to 0. The operation of UPVEC is the same except that the 
data is not read from the tape, but is already in core. Also, in UPVEC, 

NW is set equal to NEL first and then the unpacking proceeds as in the 
case of UPVECT. 

10. COMMENTS 


None. 
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LISTING 


A listing of the program is attached at the end of Section A3. 
12. TESTS 


The program has been tested with known, packed test patterns and 
also in conjunction with the routine APES discussed in Section II-A1 and 
has been found to work satisfactorily. 

\ 


102 



A3. Determination of the Part of Interest in a 'Record 


1. NAME 

Deckname: BECLMN 

Subroutine name: BECOL 

2. PURPOSE 

To find two points in a packed record such that all pixels before 
the first point and after the last point have densities equal to 63. 

3. CALLING SEQUENCE 

CALL BECOL (IX, N, II, 12, K) 

Where IX is an integer input vector consisting of N words, II, 12 
and K are output integers (see paragraph 9 below). 

4. INPUT -OUTPUT 

The input and output are through the calling sequence only. 

5. EXITS 

No abnormal exits. 

6. USAGE 

The program is written in FORTRAN IV. It is presently implemented 
on IBM 7094. 

7. EXTERNAL INTERFACES 

7.1 System subroutines: SYSLOC 

7.2 Other programs called: UNPINT 

7.3 External storage: None 

8. PERFORMANCE SPECIFICATIONS 

8.1 Storage: 260g core locations 

8.2 Execution time: TBD 

8.3 I/O load: None 

8.4 Restrictions: None 
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METHOD 


First of all, II, 12 and K are set equal to 0. Next, II is deter- 
mined so that 

II - Min { ij 3 < I < N and IX(I-) ^-377777777777 g \ 


That is, II is the index of the first word of IX, all of whose 36 
hits are not "set". If such an II cannot be found, then the program re- 
turns with II =12 = K = 0. If such an II is found, 12 is then deter- 
mined such that 


12 = Max | I j 3 < I < N~3 and IX(I) f -377777777777 g 


(Here, only I€ [3, N-3] is considered to allow for scanner label words 

Next, IX(Il) and IX(I2) are unpacked into 6 integer words each. 

Of the 6 integers resulting from IX(I1) , let the Jlth word be the first 
which is less than 63. Also, of the 6 integer words from IX(I2), let 
the J2th word be the last which is less than 63. Then, 11 and 12 are 
redefined as 


11 = (II - 1) * 6 + J1 

and 

12 - (12 - 1) * 6 + J2 

and K is defined by 


K = 12 - II + 1. 


Thus, if IA denotes a 6N vector obtained by unpacking IX, then 
IA(Il) and IA(I2) are the first and last words of IA which are less than 
63 and, hence, II and 12 mark the region of interest in the record. 

Also, K gives the length of the part of IA which is in the region of 
interest. 

10. COMMENTS 


None. 

11. LISTING 

A listing of the program is attached at the end of this section. 


104 



12 


TESTS 


The program has been tested in conjunction with the routine APES 
(Section II -Al) and has been found to work satisfactorily. 
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B. Transposition and 90° Rotations of Image Data 

V ; -***.•, ‘fj ... • • 

1. NAME 

Deck name : MTROT 

Entries: MTRANP and MROT 

2. PURPOSE 1 • 

To generate transposition or 90° rotations (clockwise or counter- 
clockwise) of an MXN data array stored on tape as M records of N words 
each (N + 1 words, in fact, the first word being the record number), when 
M <. 1024 and N < 1024. 

3. CALLING SEQUENCE 

For transposing: - > 

CALL MTRANP (NR, NC, IX, NTAPI, NIAPO) 

For 90° rotation: 

CALL MROT (NR, NC, IX, NTAPI, NIAPO, IROT) 

where 


NR a Number of rows in the input array. 

NC 53 Number of columns in the input array.. ; 

IX = Work array which can be real or integer. - 1 

NTAPI = Input tape unit. - >•' 

NIAPO = Output- tape unit. 

IROT = +1 for counterclockwise rotation and -1 for clockwise rotation. 

4. INPUT -OUTPUT 

4.1 Input: The input array should be on logical unit NTAPI. The 

records should be in FORTRAN binary format, each record containing the 
record number followed by NC words. 

4.2 Output: The output array will be written on unit NTAPI in 

FORTRAN binary format, each record continuing the record number and 
followed by NR words of data. 
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4.3 File Storage: . The number of scratch tapes used by the program 
depends on the size of the matrix to be transposed. In fact, 


n = Number of scratch tapes 85 Max j [(NR-1)/128] , [(NC-1) /128]} + 1 

* 

The first n of the tape units 2, 3, 4, 8, 10, 11, 12, 13 will be needed by 
the program. The mnemonics (B3, A4, ...) for the tape units needed will 
be printed on-line and the program will pause. The units requested should 
be provided before restarting the program. 

5. | EXITS ■ ”, . 

No non-standard exits. 

6. USAGE ' 

The program will run on any computer system with FORTRAN IV or a 
higher level of FORTRAN. It is presently implemented on IBM 7094. 

7. EXTERNAL INTERFACES 

7.1 System subroutines: SYSLOC, BSFILE 

7.2 Other programs called: RQSTWT, RWNDVJT , FLIPV, ASMBL (see 

listings attached) 

7.3 External storage used: None , . 

8. PERFORMANCE SPECIFICATIONS , 

8.1 Storage: 40770g .locations of core 

8.2 Execution time: The execution time depends on the size of the 

array to be transposed or rotated. Most of the time involved is for 
input/output operations and the time depends mainly on. the minimum number 
of 128 x 128 partitions which can cover the given array. The time taken 
for some of the arrays which were transposed and rotated using this pro- 


gram 

are shown below. 

Time 

(Sec.) 


Array Size 

Transposition 

90° Rotation 


100 x 266 . 

' 80.4 , 

. 60.9 , 


250 x 900 

296 , - t 

.. .309 


850 x 900 

1020 

1050 

9. 

METHOD „ , . : . /, 




See Section II-B of Part I. 


117 



10 


COMMENTS 


None. 

11. LISTING 

A listing of this program and the programs called by this program 
(listed in paragraph 7.2) are attached at the end of this section. The 
comments in the listings are sufficient to describe the operation of these 
subroutines. 

12. RESTRICTIONS 

The array sizes which can be handled are restricted only by the 
number of work tapes available. 
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C. Reflection of a Picture in its Horizontal Edges 


1. NAME 


Deck name: RELCT 

Subroutine name: REFLC 

2. PURPOSE 


Given input picture data on a tape, to generate a file containing 
first the reflection of a given number n of data records in the first 
record, next the data records themselves and finally the reflection of 
n records of data in the last record. 

3. CALLING SEQUENCE 

CALL REFLC (W, NR1,. NEL) 

where W is a work array which can be integer or real and NRl and NEL are 
input integers. The array W is dimensioned (NEL, NRl). (See the com- 
ments in the listing attached at the end of this section.) 

4* INPUT -OUTPUT ; 

4.1 Input: Input should be on logical unit NTAPI specified in 

the common block /INREFL/. (See paragraph 7.3 below.) 

4.2 Output: Output will be on logical unit NTAPO specified in 

the common block /INREFL/. 

4.3 File storage: The number of work tapes needed is given by 

NWTP = (NRECR - l)/NRl 

where NRECR is the number of records to be reflected in the top (and 
bottom) edge of the picture. The work tapes needed will be requested 
by an on-line message to the operator followed by a pause. 

5. EXITS 

If the number of work tapes needed exceeds 8, then the program 
will print the message :, ARRAY TOO IARGE TO HANDLE"- and return. 

6. USAGE 

This program is written in FORTRAN IV and presently implemented 
on IBM 7094. 
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7 


EXTERNAL I NT ERF AC ES 


7.1 System subroutines: SYSLOC, BSFILE 

7.2 Other programs called: RQSTWT (See listing attached at the 

end of Section II-B.) 

7.3 External storage used: The common block 

COMMON/ INREFL/NREC , NRECR , NTAP I , NTAPO 

is used to provide the inputs. NREC is the number of input records. 
NRECR is the number of records to be reflected. NTAPI and NTAPO are 
respectively the input and output logic units. 

8. PERFORMANCE SPECIFICATIONS 

8.1 Storage: 1633g locations of core 

8.2 Execution time: The execution time is highly dependent on 

NREC, NRl , NEL and NRECR. The value of NRl should be taken as large as 
possible so that the number of work tapes needed is minimized. The time 
taken was approximately 107 seconds for a test case where NREC = 1024, 
NRl = 25, NEL = 333 and NRECR = 49. For a case with NREC - 1707, 

NRl = 43, NEL * 366 and NRECR = 80, the program took approximately 164 
seconds. 

8.3 I/O load: None 

8.4 Restrictions: The size of the region which can be reflected 

(i.e., NRECR*NEL) is limited by the number of work tapes available and 
the core size. In the present setup where a maximum of 8 work tapes is 
available, NRECR*NEL will be limited to approximately 128000. 

9. METHOD 

See Section II-C of Part I. • 

10. COMMENTS - 

None. 

H. LISTING 

A listing of the program is attached at the end of this section. 
12. TESTS 


The program has been tested using a test pattern and also several 
picture data. Printouts of records and also images generated after pack- 
ing the output indicate satisfactory operation of the program. 
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D. Packing and Centering for Image Reconstruction 


1. NAME 

Deck name: PICDM 

Entries: PICDIM and PICWEM 

2. PURPOSE 

« 

To pack integer picture data (between 0 and 63) such that a speci- 
fied number of frames will be centered on a 1024 x 1024 pixel area (for 
display) in the case of entry PICDIM. In the case of entry PICWRM, the 
frames will be centered only in the horizontal direction so that the 
tape prepared can be used for writing on a film with specified width and 
writing interval. 

3. CALLING SEQUENCE 

CALL PICDIM (IB, IW, IWP) 


or 

CALL PICWRM (IB, IW, IWP) 

where IB, IW and IWP are work arrays to be dimensioned as described in 
the listing at the end of this section. 

4. INPUT -OUTPUT 

4.1 Input: Input data should be on logical unit NTPI (see para- 
graph 7.3 below). Each input record is assumed to consist of the record 
number followed by NC integer words of data where NC is defined through 
a common block (paragraph 7.3). The input records should be in FORTRAN 
binary form. 


4.2 Output: The output data will be written on logical unit 

NTPO in non-FORTRAN binary format. Under entry PICDIM, an end of file 
mark is written at the end of each file of output. Under entry PICWRM, 
no EOF is written by the program. Note that the buffer size for unit 
NTPO should be at least NDM where NDM is defined in Section 9. 

4.3 File storage: The number of work tapes needed equals NFH 

whenever NFH > 1. No work tapes are needed if NFH =1. When NFH < 1, 
the first NFH units which are not the same as NTAPI or NTAPO will be 
chosen out of the logical units 2, 3, 4, 8, 10, 11, 12, 13, 14, 15. The 
mnemonics (B3, A4, ...) corresponding to the work tapes needed will be 
printed on the line printer and the program will pause to enable the 
operator to provide the work tapes. 
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EXITS 


If the data requested to be packed exceeds the available screen 
area (under entry PICDIM) or the film width (under entry PICWRM) , then 
the program exits giving the on-line and off-line message "REQUESTED 
PICTURE EXCEEDS ALLOTTED FILM OR SCREEN AREA." If the number n of hori- 
zontal frames requested to be packed is greater than 8, the program exits 
with the message "NUMBER OF HORIZONTAL FRAMES SHD BE .LE. 8. BUT IT IS 
GIVEN TO BE n. ” 

6. USAGE 

The program is written in FORTRAN IV and presently implemented on 
IBM 7094. 

7. EXTERNAL INTERFACES 

7.1 System subroutines: SYSLOC, SKFBIN, WRITER 

7.2 Other programs called: RQSTWT, SVSCI, PCKINT, RWNDWT . The 

listings for RQSTWT and RWNDWT are attached at the end of Section B. The 
program SVSCI sets all components of an integer vector equal to a given 
integer. The routine PCKINT packs the last six bits from each of the 
words of a six -word array to form a single word. 

7.3 External storage used: The common block 

COMMON/ IP ICD/NF ILE , NR , NC , NELVER , NELHOR , NF V , NFH , 

NTPI, NTPO, IEF , ICOMP ,FILMW,WRINT,FRAMAR 

supplies the input parameters needed for the programs. 

NFILE -= Number of output files needed. 

NR = Number of rows (records) in each of the frames to be packed. 

NC = Number of columns (unpacked words per record) in each of 
the frames to be packed. ' - 

NELVER = Number of times each of the rows is to be repeated 
while packing. 

NELHOR = Number of times each of the pixels is to be repeated 
in a packed record. 

NFV = Ntanber of vertical frames to be displayed. 

NFH = Number of horizontal frames to be displayed. 

NTPI = Logical unit number of the input. 
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NTPO = Logical unit number of the output tape. 

XEF = 1 if the input frames are separated by EOF's. 

= other integer if the input frames are not separated by 

EOF’s. . , ■ . . ' : 

ICOMP = 0 if the true (positive) image is to be packed. 

■ = 1 if the complement (negative) image is to be packed. 

FILMW = Width of film in millimeters. 

WRINT = Writing interval in microns. 

FRAMA.R = Vertical margins needed between frames in millimeters. 

The last three parameters need be used only for PICWRM. 

Note: The number of records and the number of words per 

record in each of the frames to be packed during one 
call should be the same for all the NFILE*NFH*NFV 
frames of data. 

8. performance SPECIFICATIONS 

8.1 Storage: 1163g locations of core.. ' > 

8.2 Execution time: The execution time depends on the sizes of 

picture data arrays to be packed. If NFH (number of frames in the hori- 
zontal direction) equals 1, the program involves only packing of data 
and writing it on the output tape. If NFH > 1 then the data will be 
written on work tapes, the work tapes are rewound and then the output 
records are written. Thus, for a given product NFH*NFV, NFH = 1 results 
in faster implementation than NFH > 1. (Of course the desired value of 
NFH depends on the user's needs.) Two typical times for the case NFH = 1 
are shown below: < ~ . 

NFILE = NELVER = NELHOR = NFV = NFH = 1 in both cases. I 

For NR = 900, NC = 850, PICDIM was executed in approximately 
170 seconds. 

For NR = 1707, NC = 366, FILMW = 25.0, WRINT = 12.5 and FRAMA.R = 0.0, 
PICWRM was executed in approximately 130 seconds. 

8.3 I/O load: None 

8.4 Restrictions: NFH <. 8 (due to limit on number of available 

work tapes) 
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METHOD 


9.1 Entry PICDIM: The number of packed words needed to produce 

one line of display is 171, corresponding to 1026 unpacked words. Using 
this fact, the horizontal margin available is computed in terms of num- 
ber of packed words. The margin is equally distributed between frames 
and the two ends. Similarly the vertical margin available is computed 
as (1024 - NR*NELVER*NFV) and distributed equally between frames and the 
two ends. Denote by* MARGVE the vertical margin between the frames. 
Denote by MARGHO the horizontal margin between the frames. 

The first NFV frames of the picture are packed and written on the 
first work tape, the second NFV frames on the second tape, and so on. 

Each record on the Ith work tape will have the Jth word equal to 
zero for all J such that 

i ' 

J 4 [ I * MARGHO + (1-1) + 1, I * MARGHO + NHOr] 

where 

NHOR = (NC * NELHOR) /6 


After writing the packed records corresponding to NFV * NFH frames 
on the NFH work tapes, the work tapes are rewound. The Ith record from 
each of the ; work tapes is read and the NFH such records are added to 
form the Ith record of output, for 1=1, 2, 3, ... . 

I * 

. If ICOMP = 1, the data values are subtracted by 63 before packing. 

* ' ' 

: .9.2 Entry PICWRM: In this case the procedure is the same as in 

entry PICDIM except that 

MARGVE = 1000 . *FRAMAR/WRINT 

and the number of" packed words per output record is given by 
NDM = 1000.*FILMW/(WRITN*6.) + 1 

and the total horizontal margin is computed using NDM instead of 171. 

10. COMMENTS 


The margin areas appear light (0 density) regardless of whether 
ICOMP = 0 or 1. 
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LISTING 


A listing of the program is attached at the end of this section. 
12. TESTS 


The program has been tested for a large number of cases including 
test patterns and found to work satisfactorily. The pictures shown in 
this report were all prepared using this program. 
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SECTION III. IMAGE PROCESSING PROGRAMS 


Al. Filter Design by Helms' 4-T*s Method 


1. NAME 

Deck name : HELM 

Subroutine name: HELMS 

2. PURPOSE 

This subroutine generates a finite number of filter weights which 
approximately yield a specified frequency response function. Helms' 
4-T's method is used (hence the name). 

3. CALLING SEQUENCE 

CALL HELMS (IF LAG, N, SAFR, TIMR,L,LMAX, ERROR) 

where 


IFIAG ° An integer between 1 and 6 indicating the type of filter. 
(1: Low Pass; 2: High Pass; 3: Band Pass; 4: Low Emphasis; 5: High Em- 
phasis; 6: Band Emphasis) 

N = Number of samples of the frequency response. 

SAFR = Complex output array dimensioned N which will contain the 
frequency response samples. 

TIMR ® Complex array dimensioned N which handles the impulse 
response (i.e., filter weights) corresponding to SAFR. 

L ° Input constant which gives the initial value of the number of 
filter weights with which the approximation to desired frequency response 
is to be attempted. Also, output giving the final number of filter 
weights used for the approximation. 

IMAX = Input constant which gives the upper limit on the number of 
filter weights. 

ERROR = Permis sable error “ MAX | Gj - Gi I ; (an input 

i-l,...,N 1 

constant) where G^ and G| are the ith components of the given and the 
approximating frequency responses, respectively. 
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INPUT-OUTPUT 


4.1 Input: Input to HELMS is through the calling sequence only. 

4.2 Output: L, ISMALL and the first L components of the real part 
of TIMR are written on logical unit 11. 

L = Number of filter weights (integer) . 

ISMAIL = An integer representing the shift in the approximating 
impulse response (see Section 9). 

TIMR = Complex array approximating impulse response. 

L, ISMALL are written as one record and Real parts of the components 
of TIMR are written as L records of one word each. 

For each attempted value of L, the value of L and ERR, the error 
in approximation are printed off-line so that the user knows whether the 
program stopped with L = LMAX and/or with ERR < ERROR. 

4.3 File storage: None 


5. EXITS 

The program exits when ERR < ERROR or L >_ LMAX during the itera 
tions (see Section 9). If the exit occurs due to the latter condition, 
L will be set equal to LMAX and the corresponding impulse response will 
be written on unit 11. 

6. USAGE 


This program will run on any computer system with FORTRAN IV or a 
higher level of FORTRAN. It is presently implemented on IBM 7094. 

7. EXTERNAL INTERFACES 

7.1 System subroutines: SYSLOC 

7.2 Other programs called: SAMPLE, TRIM (Sections A3 and A4) and 

FFOURT (a Fast Fourier Transform routine) 

7.3 External storage used: Common blocks 

/WOR/WORK(1024) ,W0RK1 (1024) 

/ PARMTR/ OMEGC 1 , OMEGT 1 , 0MEGC2 , 0MEGT2 , P 1 , P2 , A 

Common block /WOR/ is for complex work arrays which can be made 
common with the main program (see Section 10). 
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OMEGCl,OMEGTl,OMEGC2,OMEGT2,Pl,P2,A are constant parameters defining 
the frequency response to be approximated. ; r 

8. PERFORMANCE SPECIFICATIONS 

8.1 Storage: 10505g locations of core. 

8.2 Execution time: The execution time depends on the number of 

samples (N) and the number of trials required to obtain the desired 
approximation to the transfer function. Since the filter lengths are 
doubled at each successive trial and the length is bounded by LMAX, the 
number of trials is less than or equal to log 2 (LMAX/L) where L is the 
initial filter length specified. It has been found that the average 
time per trial is 169 seconds where N = 512 and it is approximately equal 
to 43 seconds when N = 256. 

8.2 I/O load: The maximum number of lines of (print) output due 

to HELMS is the smallest integer greater than log 2 (LMAX/L) where L refers 
to the input value of L. 

8.4 Restrictions: N < 1024 and (Input L) LMAX <1024 

Dimensions of SAFR and TIMR in the calling program should be 
greater than or equal to N. 

(Even though the program works for all LMAX < 1024, it is advis- 
able to use LMAX < 310 so that the filter obtained can be used in con- 
junction with FFTCNV.) 

9. METHOD 

Helms' 4-T's method (Transform-Truncate^Transform-Test) is used for 
approximating the assumed frequency response., 

a. The assumed frequency response is sampled and N equally spaced 
(complex) samples G n are generated as described under paragraph 9 of the 
documentation of subroutine SAMPLE (Section A3). For the purposes of 
filtering,' it may be assumed that the weights desired are real numbers 
and hence the frequency response assumed must have complex conjugate sym- 
metry. That is, G(-w) = G*(w) for all w. 

b. The inverse DFT {g n | of the sequence |G n | is found. 

c. Ah integer i is chosen such that 
0 < i < N-l 

and m 

E ' E I S k | 

k=i 
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is the minimum where 

m = (i+N-L-1) (Mod N) 


d. A new sequence jg^j is defined as follows 


and 


3 n 


0 for n = i, i+1 i+N-L-1 (Mod N) 


Re(g n ) otherwise 


The real parts of g n are taken in the above equation since the filter 
weights needed are real and the sequence jg n } has, in general, nonzero 
imaginary parts due to computational errors. 

e. The DFT {g^} of jg'j is found and, the error in approxi- 
mating is computed, where 


E 1 = 


Max |G n > g' 

0 < n <. N-l 


It is easy to see from the definition of DFT that 

. N-l 


»\ T.jnk 


Max (g k - g£) W 1 

0 <. n <. N-l 

k=0 


m 


Max |y^ g R W nk + 

0 <. n < N-l 1 7“^ 

- k=i 


k $ [i,m] 


nk 


“a* \y M i»" k l + y |i»(s k )| |w nk 

0 < n < ^ 


where 


W = exp (-2 j/N) 
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Therefore , 


E i ^ 




m 


N-l 

Max 

N-l 

£ 

| S k| 

+ £ 

n _< 

-k=i 


k=0 




i « j 


1 — E + ^ [ ImC8 k ) | 


k=Q 


Hence, the minimization of E defined in step (a) implies minimization 
of an upper bound on the maximum error in approximating the frequency 
response. 


f. Having computed E-^, it is checked whether E-^ < ERROR (the 
pre-specified permissible error) or L UMAX. If neither of the condi- 
tions holds, then the value of L is doubled and the steps starting from 
(c) are repeated. If L becomes greater than LMAX when L is doubled, 
then L is set equal to LMAX and the steps starting from (c) are repeated. 

g. Thus, at the conclusion of step f either E^ < ERROR or 

L = LMAX (or both). Next, a new sequence |g^| is defined by shifting 
all the (N-L) contiguous zero values of the sequence | to the 

last (N-L) positions. That is, 8 


II 

0 

= 

g i+N 

It 

1 

• 

g I+N 

.11 

L-2 

• 

e 

g i-2 

it 

L-l 

= 

8 i-l 

.it 

L 

= 

0 > 

it 

L+l 

• 

0 

* 

ii 

N-l 

• 

• 

S3 

0 > 


> 


L terms 


N-L terms 
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where the subscripts are taken Modulo N. Now, g«, g^_-^ are the 

filter weights of a "Finite Impulse Response (FIR)" filter whose fre- 
quency response approximates the given frequency response except for a 
linear phase difference which is completely specified by specifying the 
number i denoted by ISMALL in the program. Since i indicates the shift 
in the filter weights (actually, the weights g^ are shifted forward by 
L-i subscripts to give weights gjj) , the number i is called a "shift 
parameter". 

h. The values of L, i and gQ, ..., g" are written on tape as 
the outputs of the subroutine HELMS. 

A flow chart is shown in Figure 3.1 showing the various steps in 
the subroutine. 

10. COMMENTS 

At the end of a call of HELMS, the work arrays WORK, and W0RK1 will 
contain the approximating frequency response and the exact impulse 
response (i.e., the IDFT of the Sampled Frequency Response) , respectively. 
By including a common block /WOR/ in the calling program WORK and W0RK1 
can be retrieved if necessary (e.g., to get plots of original and approxi- 
mating frequency response). 

11. LISTING 

A listing of the program is attached at the end of this section. 

12. TESTS 

The program has been used for generating several low pass, high 
pass, band pass and low emphasis, high emphasis and band emphasis filters. 
See Section III.C of Part I for some examples of the original and approxi- 
mating filter transfer functions. 
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Figure 3.1: Flow Chart for Subroutine HELMS (Continued) 
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A2. 


"Optimum” Filter Design by Helms' 4-T's Method 


1„ NAME 

Deck name: HELM2 

Subroutine name: HELMS2 - . 

2. PURPOSE 

This subroutine generates a filter with the smallest number of filter 
weights which approximate a specified frequency response number to within 
a specified error. 

3. CALLING SEQUENCE 

CALL HELMS2 (IFLAG, N, SAFR, TIMR, L, LMAX, ERROR) 

where all the calling arguments except L have the same meaning as in the 
case of HELMS (Section Al) . The argument L is only an output giving the 
final value of the number of filter weights. 

4. INPUT -OUTPUT ■ ■■ 

Same as in the case of HELMS. 

5. EXITS 

No abnormal exits. ; 

6. USAGE 


The program is written in FORTRAN IV. It is presently implemented 
on IBM 7094. 

7. EXTERNAL INTERFACES 
Same as for HELMS. 

8. PERFORMANCE SPECIFICATIONS 

8.1 Storage: . 1356 q’ locations of core ■* • 

* i 

8.2 Execution time: It can be seen from paragraph 9 that the num- 

ber of trials of filter lengths is less than or equal to log 2 (LMAX - 2) + 3. 
The time per trial is the same as for HELMS. 

8.3 I/O load: Same as for HELMS. 

8.4 Restrictions: Same as for HELMS. 



9. METHOD 


This program differs from HELMS only in the way a new L is defined 
in step f (paragraph 9, Section Al) . Denote by E^ the values of 
filter length and error, respectively, during the ith trial. The program 
tries = 1 and L 2 * LMAX. It is assumed that Ei > ERROR. (It is 
reasonable to assume this because most filters of interest cannot be 
approximated with just one filter weight. One filter weight is just a 
constant gain.) If E 2 > ERROR, the program writes LMAX filter weights 
on unit 11 and returns. If E 2 ERROR, the following rule is used for 
defining L i+ ^ from L^ , L 2 , ...» L^ 

L i+1 = (Lj+L.m ■ 

where 

j = Max ) k < i | (Ei - ERROR) (E k - ERROR) <_ 0 } 

| ■■■!■ ■ ■ ■■ ■ : ' r . ;• 

The iterations are stopped when I^j_i and I^.have been found such that 
1^ " ^m-l = 1* Then, the final value of L will be chosen from the 
sequence Li, . . „ such that the smallest value of Li with Ei <. ERROR 
is obtained. 

It is quite simple to verify that the number m of iterations needed 
is limited by 

m < log 2 j LMAX - 2 j + 3 « 

iO. COMMENTS 

: Same: as for HELMS. r- [ ,i 


11. LISTING 

A listing of the program is attached at the end of this section. 

12. TESTS 

This program was tested by designing several band pass and band 
emphasis filters. With N = 256, the number of iterations was found to 
be equal to 10. 
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IF(ER(I).GE.ERRJR.AND.ER(I).LE.ERMN) ERMNl=ER(I) 

— T F ( ER ( I )„GE. ERROR.iND. f ER (I) ,TT. ER.MN.OR, ER (T)„ EQ. ERMN. SNDTl'l IT) .TT 
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END 



A3. Sampled Frequency Response Generation 


1. NAME 

Deck name: SAMPL ' 

, Subroutine name: SAMPLE 

2. PURPOSE 

To generate a given number of equally spaced samples of a specified 
frequency response function. 

3. CALLING SEQUENCE 

CALL SAMPLE (IFLAG, N, SAFR) 

where 

IFLAG = An integer input between 1 and 6 indicating the type of 
filter. (1: Low Pass; 2: High Pass; 3: Band Pass; 4: Low Qnphasis; 

5: High Qnphasis; 6: Band Emphasis) 

N = Number of samples of the frequency response (input integer). 

SAFR = Complex output array dimensioned N which will contain the 
frequency response samples. 

4. INPUT -OUTPUT 

4.1 Input: Input to SAMPLE is through the calling sequence only. 

4.2 Output: Output is through the calling sequence only. 

4.3 File storage: None. 


5. EXITS 

No non-standard exits. 

6 . USAGE 


This program will run on any computer system with FORTRAN IV or a 
higher level of FORTRAN. It is presently implemented on IBM 7094. 

7. EXTERNAL INTERFACES 

7.1 System Subroutines: SYSLOC 
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7.2 Other programs called: FREDMl, FREDM2, FREDM3, FREDA (See 

listings at the end of this section) 

7.3 External storage used: Common block 

COMMON/PARMTR/OMEGC 1 , 0MEGC2 , 0MEGT2 , PI , P2 ,A 

contains parameters describing the frequency response (to be defined in 
the calling program). 

8. PERFORMANCE SPECIFICATIONS 

8.1 Storage: 170g locations of core 

8.2 Execution time: TBD 

8.3 I/O load: No input /output commands. 

8.4 Restrictions: None 

9. METHOD 


The expressions for the magnitude of the given band limited frequency 
response are assumed to be defined by the functions FREDMl, FREDM2, FREDM3 
and the phase of the frequency response is assumed to be defined by func- 
tion FREDA. (FRED = Frequency REsponse Definition; M s Magnitude; 

A = Angle) 

In the present implementation, FREDMl defines a low pass character- 
istic |G^(u)| shown in Figure 3. 2a. In the interval [- u s /2, w s /2 ] , 


G x (CJ) 


1 for 1 0 > | < 

0 for { wj > OJ-j 


and 


where 


( Utrj, _ j(*>| \P 

[(Jy -<J C j f0r ^C < | ^ W T 


"s = 

f s = 

Ax 

U) c , <»> T 


sampling frequency = 1/Ax 
sampling interval 
and p are specified parameters 
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Function FREDM2 defines a high pass characteristic G 2 (u>) shown 
in Figure 3.2b. 


G 2 (CJ) = 1 - G^(CJ) for 

Function FREDM3 defines a band pass characteristic G^Cw) shown 
in Figure 3i2c. 


G 2 ((x>) - G^ 2 (oj) ~ 

where 


G li(tL») = G^Cd) with (j c = o> ci , (J T = 0» Ti and p » p £ 
for i = 1, 2. 

It is assumed that the calling program specifies the parameters w , 
w as fractions of o • It is then to be noted that i 

1 s 

0 < 0MEGC1 < 0MEGT1 < 0MEGC2 < 0MEGT2 < 0.5 


in the case of the band pass filter. Further, OMEGC1 and OMEGTl are 
specified as 0 in the case of low and high pass filters. 

Now, sample numbers n = 1, ..., N are assigned to frequencies between 

- <j /2 and u / 2 as follows: 
s s 

= [(n-1) /N ] Cd g /2 for n = 1, ...» N/2 

and 

W n = - [(N - n+1 ) /n] Djl for n = N/2 + 1, ..., N 

The SAmpled Frequency Response, SAFR, is defined by 

SAFR(n) = | G £ ( <4> n ) j exp(j0(<d n ) ) 

where 0(w) is the phase defined by function FREDA and i is chosen to 
be 1, 2, 3 for low, high, band pass, respectively. 
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If IF LAG is specified as. 4, 5 or 6, then low, high, or band emphasis 
filter characteristics are obtained by modifying the corresponding "pass" 
characteristics. Only the magnitude characteristics are changed such that 
the "stopband" gain is A instead of zero and the "passband" gain is 1. 

In f ac t , 


' SAFR(n) = [ | G i(4d n )| (1-A) + Aj exp(j 0(6> n ) ) 

where i = IFLAG-3 and the other symbols are as defined above. 

10. COMMENTS 

The polynomial approximation to the low pass characteristic G^(w) 
can be replaced by any other desired type of approximation by modifying 
the function FREDMl suitably. 

11. LISTING 

A listing of subroutine SAMPLE and function subprograms FREDMl, 

FREDM2 , FREDM3 and FREDA are shown at the end of this section. 

12. TESTS 

The program has been tested in conjunction with HELMS. See paragraph 
12 under Section Al : . 
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A4. Impulse Response Truncation 

1. NAME 

Deck name : TRUIM 

Subroutine name: TRIM 

2. PURPOSE 

To truncatei the" impulse response (or any complex vector with N com- 
ponents) by setting (N - L) contiguous components to zero, considering 
the impulse response to be periodic with period N. The components to be 
set equal to zero are chosen so that the sum of their magnitudes is mini- 
mized. 

3. CALLING SEQUENCE 

CALL TRIM (TIMR, N, L, ISMALL) 

TIMR “ Input vector with N complex components. Also output vector 
whose (N - L) components are equal to the corresponding components of" the 
input vector. 

N = Input integer 

L = Input integer equal to the number of components of TIMR to be 
retained. 

ISMALL “ Output integer equal to the subscript starting from which 
(N - L) components of TIMR are set to zero. 

4. INPUT-OUTPUT 

4.1 Input: Input is through the calling sequence only. 

4.2 Output: Output is through the calling sequence only. 

4.3 File storage: None. 

5. EXITS 

No non-standard exits. 

6. USAGE 

This subroutine will run on any computer system with FORTRAN IV or a 
higher level of FORTRAN. It is presently implemented on IBM 7094. 
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7 


EXTERNAL INTERFACES 


7.1 System subroutines: SYSLOC 

7.2 Other programs called: None 

7.3 External storage used: None 

8. PERFORMANCE SPECIFICATIONS 

8.1 Storage: 335g core locations 

8.2 Execution time: TBD 

8.3 I/O load: No I/O commands 

8.4 Restrictions: None j 

9. METHOD 

Let Xq, . . „ , x N ^ be the components of the input vector. Then,’ 
for k = 0, N-l , the program computes 

■wf 


L-l 



where 


X k+i 


x with m = (k+/ ) (mod N) 
m 



The value of k for which s^ is minimum is determined. The components 
x i> x i+l» •••» X i+N-L are set t0 zero » ^ he subscripts being again taken 
modulo N. 

10. COMMENTS 
None. 

11. LISTING 

A listing of the program is attached at the end of this section. 

12. TESTS 


The program has been tested in conjunction with HELMS (see para- 
graph 12 under Section Al. 
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Bl. Point Operations - General 


■1. NAME 
SCALE 

2 . PURPOSE 


To scale a set of numbers linearly or nonlinearly so that the data 
lie between 0 and 63 to make the data suitable for display. The program 
handles real as well as integer inputs, and as many data files as speci- 
fied. 

3. CALLING SEQUENCE 

CALL SCALE (NFILE, NR, NC, IW, RW, PAR, FUNC, IH) 


where 


NFILE = Number of input files 
NR = Integer input array where 

file. 


(Integer <. 20) 

NR(I) = number of records in the Ith 


NC = Integer input array where NC(I) = number of words of picture 
data per record of Ith file. (It is assumed that each record of Ith file 
consists of NC(I) + 1 words, the first word being the record number.) 

IW and RW are work arrays for handling the input and output records. 
(See comments in the attached listing.) 

PAR is a real array consisting of parameters defining the external 
scaling function FUNC. 

IH is an integer array for storing the histogram of the output data. 

4. INPUT-OUTPUT 

See paragraph 7 for the definition of variable names used below. 

4.1 Input: Input data are on unit NTAPI. The first word in each 
record is an integer (the record number) and the remaining NEL words are 
integer or real data. (If real, ITYPE should be specified as 1.) It is 
assumed by this subroutine that the input data are in FORTRAN binary for- 
mat. 


4.2 Output: The output data appear on unit NTAPO as integers. 

The record number is written in the beginning of every record and NEL 
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data words follow. The format is FORTRAN binary. If input files are 
separated by EOF (i.e., IEOF is given as 1), the output files also will 
be separated by EOF. 

4,3 File storage: None 


5 . EXITS 


No non-standard exits. 

6. USAGE 

This program' will run on any computer system with FORTRAN IV or a 
higher level of FORTRAN, It is presently implemented on IBM 7094. 

7. EXTERNAL INTERFACES 

, \ 

7.1 System subroutines: SYSLOC, SKFBIN 

7.2 Other programs called: FUNC (an EXTERNAL function to be speci- 

fied by the user in the form FUNC (PAR, X). See Section 3, above.) and 
IRHSTG (see Section B2, below). 

7.3 External storage used: Three common blocks are used: 

COMMON/lNSCAL/NTAPI, NTAPO , ITYPE, IEOF , IMN, IF LG , IHST 
COMMON/FILE/IF ILE 
COMMON/WSCAL/RWB , RWS , RWM 

The block /INSCAL/ supplies the input parameters to the program. 
The meanings of all the variable names are given in the listing attached. 
The block /FILE/ is useful if the scaling function FUNC is to be made 
independent on IFILE, which is a do loop index giving the file number 
being scaled. In that case, the block /FILE/ should be made common to 
SCALE and FUNC. The block /WSCAL/ consists of work arrays RWB, RWS, RWM 
which store the maxima j minima and means of the files being scaled (under 
the option IFLG f 1 and IMN f 0) . This block should be made common to 
SCALE and the calling program if the maxima, minima and/or mean are to be 
used again. A typical example of the use of this common block is when 
it is necessary to scale the same file with several scaling functions, all 
of which depend upon the maximum and minimum of the input data. In this 
case, one would call SCALE first with IFLG f 0 and IFLG ^ 1 (say IFLG = 2) 
and perform scaling with the first function and for all the successive 

functions one would call SCALE with IFLG = 1. 

* , 

8. PERFORMANCE SPECIFICATIONS 

8.1 Storage: 665g locations of core. 

8.2 Execution time: The execution time is highly dependent on NFILE, 

NR, NC, the option used (see comments in listing), and the scaling function 
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FUNC. For a case where NFILE = 1, the array size was 1707 x 366, linear 
scaling' was used and the option used was such, that the minimum and the 
maximum of input data were found before scaling, it was found that the 
time taken was 400 seconds. For the same data when linear scaling was 
used with threshold and the minimum and maximum were assumed to be given, 
the time taken was 257 seconds. 

8.3 I/O load: No printouts on-line or off-line. 

8.4 Restrictions: NFILE <. 20. 

9. METHOD 

This program has several options as detailed in the comments in the 
listing attached. The method given below applies to the case when 
IFLG f (0 or 1), IMN f 0, and IHST ± 0. 

(i) First, all the data records are read and the maximum, mini- 
mum, and mean values of data in each data file are determined 
separately and stored in work arrays RWB, RWS, and RWM, 
respectively. (RWB, RWS, and RWM are dimensioned 20 which 
permits handling of 20 files.) 

(ii) The input tape is rewound. 

(iii) Each of the files is scaled separately as follows: 

When the Ith file is being scaled 

PAR(l) = RWB(I) 

PAR (2) = RWS (I) 

PAR (3) = RWM(I) 

and the output Y corresponding to input X is given by 
Y = FUNC (PAR, X) 

The output Y is computed for each data element and rounded 
off to the nearest integer. If the integer is greater than 
63, it is set to 63. If the integer is less than zero, it is 
set to zero. The number of occurrences of the integer I in a 
given output record is added to IH(I) as soon as the record 
has been computed. 

10. COMMENTS 


The function FUNC is an arbitrary point operation. It can be chosen 
depending on the particular needs of the user. This approach has been 
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taken because "the true representation 1 * of picture data is yet an unde- 
fined phenomenon. If one has to represent data points which are known to 
be all of the same sign, it would be advantageous to use a scaling func- 
tion which is different from the one in the case where the data can take 
on both positive and negative values. Also, if the data have a very wide 
range and parts of them are to be displayed with increased contrast and 
other parts with decreased contrast, a suitable nonlinear function may be 
defined. 

11. LISTING 

A listing of the program is attached at the end of this section. 

12. TESTS 


The program has been tested for a large number of cases and the 
scaled pictures have been displayed. The filtered pictures shown in this 
report were all first scaled using this program and then packed using 
PICDIM or PICWRM (Section II-D) . 
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RETURN 

EnlU 



B2. Point Operations - Non-Negative Integer Input 


A set of simple and fast routines using a table look-up method Ha's 
been developed for performing point operations on non-negative integer 
input data. The listings of these routines are attached at the end of 
this section. Their names and purposes are tabulated below. 


Name 

| 

Other Routines 

Deck 

Subroutine 

Purpose 

Called 

TBLKP1 

. 

TBLKPV 

1 

Given a vector IX and a "table 
array" ITBL, to generate a vector 
IY such that IY(I) = ITBL(IX(I)+ 
1). 

None 

TBLKP2 

TBLKPM 

Given an input array consisting 
of non-negative integers on tape 
unit NTAPI, and a "table array" 
ITBL, to generate an output array 
on tape unit NTAP) , each of whose 
records is obtained by using 
TBLKPV on the corresponding input 
record. 

TBLKPV 

IRHIST 

IRHSTG 

Given integer arrays IX and IH 
to add the histogram of array IX 
to the array IH. 

None 

i 

IP HIST 

IPHSTG 

1 . - . - . 

To find the histogram of a 
given integer data array con- 
sisting of several records. 

IRHIST 

DISLN 

i 

DISLIN 

Given the histogram IH, to set 
up a "table array" ITBL which 
can be used in routines TBLKPV 
and TBLKPM using a particular 
type of density transformation 
which linearizes the distribu- 
tion function. 

None 

i 


The routine DISLIN is a particular transformation which serves to 
stretch contrast in some pictures in which the density range is not fully 
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utilized. The above programs can be used with any other density trans- 
formation which can be supplied in the form of a "table array" ITBL. 
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C. Non-Recursive Filters - Frequency Domain 
Implementation Using Fast Convolution 


1. NAME 


Deck name: FILT 

Subroutine name : FILTER 

2. PURPOSE 


To filter a given number of data records, each with a given number 
of words (e.g. , picture data), using a filter whose weights are given. 

3. CALLING SEQUENCE 

CALL FILTER (NREC , NEL, IFLAG, IX, IY) 


where 


NREC = Number of records (input). 

NEL = Number of words in each record (input) . 

IFLAG = Input integer equal to 1, 2 or 3 depending on how the ends 
of a, record are to be "padded" to avoid "edge effects." 

IX = Integer array into which the subroutine reads the data to be 
filtered. At the end of call, IX will have the last record of input data. 

IY = REAL array which handles the filtered data. At the end of call, 
IY will have the last record of filtered data. 

4. INPUT -OUTPUT 

4.1 Input: 

(i) Unit 11: The number of weights (L) and the shift parameter 

(ISMALL) are read as one binary record. The filter weights 
are L real numbers which are read as L binary records. 

(ii) Unit 8: IREC, (IX(I), 1=1, NEL) are read as one binary 

record, where IREC is the record number and IX is the 
IRECth integer data array. The reading takes place NREC 
times corresponding to IREC = 1, ...» NREC. 

4.2 Output: 

Unit 13: IREC, (IY(I), 1=1, NEL) are written as one binary 

record where IREC is the record number and IY is the IRECth 


# 
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- /filtered, data array. The writing takes place,»NREC . times 

corresponding to IREC * 1, NREC.. . « 

4.3 File storage: None. L 

> i ■ 

5. EXITS 

No non-standard exits. 

6. USAGE 

This program will run on any computer system with FORTRAN IV or a 
higher level of .FORTRAN. It is presently implemented on IBM, 7094. 

7. EXTERNAL INTERFACES 

7.1 System subroutines: SYSLOC 

7.2 Other programs called: FFTCNV (a fast convolution routine which 

uses Helms' "select-save" method discussed in Section III-B1 of Part I). 

8. PERFORMANCE SPECIFICATIONS 


8.1 Storage: 5233g core locations 

8.2 Execution time: The execution time depends on the filter size 

and the data array size. Typical times are listed in Table 3.5 of Part I. 

8.3 I/O load: See INPUT-OUTPUT (Section 4). 

8.4 Restrictions: The number of filter weights must be less than 

or equal to 310. 

Number of elements (NEL) per record is restricted to 2644, 
limited by the core capacity. NEL < 2141 will result in more efficient 
input/output buffer allocation. 

9. METHOD 

The routine uses a more basic subroutine called FFTCNV which imple- 
ments Helms' select-saving method of evaluating convolutions. See 
Section III-B1 of Part I for the details of array handling. 

If { x^ } , k = 0, . . . , N-l defines an input record, then the pur- 
pose of the parameter IFLAG in this program can be explained as follows: 

If IFLAG = 1, the sequence | x r| is considered to be periodic with 
period N. That is, x_-^ = x jy “ Xq, and so on. 



If IFLAG a 2, the sequence | is assumed to be symmetric about 
the end points. That is, x_^ m x^, Xjj ® x N-2» an< ^ so on * 

If IFLAG = 3, the data is assumed to be zero outside the given 
record. 

10. COMMENTS 
None. 

11. LISTING 

i . . • * • • ' • 

A listing of the subroutine is shown at the end of this section. 

12. TESTS 

See Section III-Bl for examples of filtered test patterns. 
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Recursive Filter Implementation 


1. NAME 


Deck name: TDFLl 

Subroutine name: TDFIL1 

2. PURPOSE 

To implement a "product-type” separable recursive filter which is a 
low pass filter in the vertical direction and any filter in the horizontal 
direction. 

3. CALLING SEQUENCE 

CALL TDFILl (NEL, HORFIL, PAR, IX, X, W, Y, Wl) 

where NEL is the number of words of picture data per input (and output) 
record and all other arguments are described in the comments in the at- 
tached listing. 

5 . INPOT-OUTPUT 


4.1 Input: The data to be filtered should be on logical unit 8 and 

the data including its reflections in the top and bottom of the picture 
generated by subroutine REFLC (see Section II-C) should be on logical 
unit 10. The program assumes that the number of words in each record on 
units 8 and 10 is NEL + 1, the first word being the record number and 
the remaining NEL being the data words. The input data may be either 
real or integer and should be in FORTRAN binary format. 

4.2 Output: The output data will be real. It will be written on 

unit 11 in the same format as the input. 

4.3 File storage: None 

5. EXITS 


No non-standard exits. 

6. USAGE 

This program has been written in FORTRAN IV. It is presently 
implemented on IBM 7094. 

7. EXTERNAL INTERFACES 


7.1 System subroutines: SYSLOC, BSFILE, SQRT 
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7.2 Other programs called: VADDI, VSADDI, VADDR, VSADDR, VSUBI, 

VSSUBI, VSUBR, VSSUBR (all entries under deck name V^DSUB; see attached 
listing) RLPFLl (a recursive low pass filter routine in one dimension; 
see attached listing). 

7.3 External storage used: Two common blocks are used: 

COMMON/lNPAR/NREC , NRECR, ITYPE 

COMMON/ VERPAR/LR , NORM 

See comments under listing attached for definitions of vari- 
ables in the common block. 

8. PERFORMANCE SPECIFICATIONS 

8.1 Storage: 1430g locations of core 

8.2 Execution time: The execution time depends on the size of the 

array to be filtered, the size of the filter (parameters PAR(l) and LR) , 
the horizontal filter routine HORFIL and whether normalization is re- 
quested or not. For example, a 1000 x 320 picture array was filtered 
using a recursive band pass filter routine for HORFIL and LR = 20, 

PAR(l) = 25 and NORM = .FALSE, in approximately 8.5 minutes. With 
NORM = .TRUE., the same filtering operation took about 14.2 minutes. 

8.3 I/O load; None 

8.4 Restrictions: None 

9. METHOD 

The implementation when NORM = .FALSE, follows that given in Table 

3.4 of Part I for "Product" type separable recursive filters which are low 
pass in the vertical direction. Note that L-^ and L 2 in Table 3.4 corres- 
pond to PAR(l) - 1. and LR - 1, respectively, in the program. Also, the 
division by the constants 2L^ + 1 and 2L£ + 1 in the horizontal and verti- 
cal low pass filter realizations shown in Figures 3.22 and 3.26 of Part I 
is not performed since most filtered outputs are linearly rescaled before 
displaying. When NORM = .TRUE., the output is normalized so that only 
the variations of the input with respect to the neighbors in a (2LR-1)* 
(2PAR(1)-1) is significant. The implementation of normalization in this 
case follows that shown in Figure 3.28 (Part I). To get the matched 
filter output given by Equation (3.53) of Part I, one needs to divide 
the output of this program by a constant. The value of this constant 
depends on the subroutine HORFIL. If HORFIL is either RHPFL1 or RBPFL1 
(see attached listings), then this constant is 
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(2LR-1) 



4 / 


where g^ are the corresponding filter weights. 

10. COMMENTS 
None. 

11. LISTING 


A listing of this subroutine and the subroutines referred to in this 
section is attached at the end of this section. 

12. TESTS 

The Figures 3.29, 3.30, 3.31 of Part I show the results of using 
the routine TDFIL1 with an external routine RBPFL1 in place of HORFIL. 

This is equivalent to a matched filter which enhances vertical lines. A 
typical record was taken in all these cases and normalized by dividing 

it by (2LR-1) £ gj^? . In the noise-free case, it was found that the 

outputs of the program were exactly equal to the theoretically calculated 
values. When the template (filter) matched the line in the test pattern 
with the correct width (13 pixels) the output was unity. Figures 3.3a and b 
show a record of the noise-free test pattern and a record of the filtered 
noise-free test pattern. Figure 3.3c is constructed as follows. The 
locations of the center lines of the wide vertical lines in the test pat- 
tern are determined. The filtered output at these points (i.e., when the 
center line of the template coincides with the center line of the wide 
vertical lines) is plotted against the line width. Figures 3.4a, b and c 
show the corresponding plots for the case of a noisy test pattern where 
the noise is approximately Gaussian with zero mean and a variance of 8. 

This gives a signal-to-noise ratio of .687 for the line width 13. The 
corresponding theoretically expected value of output is 

1/ yi + 1/ .687 = .636 . It is seen that this is very close to the 

value corresponding to line width 13 in Figure 3.4c. 
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CORR. 


^ONE RECORD OF NORMALIZED CORRELATIONS 



Figure 3.3b One Record of Matched Filtered 
Noise-Free Test Pattern 
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CORR. 



Figure 3.3c Filter Output at Centers of Vertical 
Lines (Noise-Free Test Pattern) 
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